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ABSTRACT 

The three-dimensional local induction problem is 
considered for various source configurations. The 
electromagnetic field solutions for layered conductors are 
developed and solutions for laterally inhomogeneous 
conductivity structures are found. A comparison of 
theoretical and analogue solutions is made. The effect of 
the source on perturbation and induction arrows is 
investigated. 

Three source field and conductivity structure 
combinations are considered. The first model of a 
uniform source above a resistive anomaly approximates an 
island in the deep ocean for Geomagnetic Bay type 
disturbances. The second model consists of an L-shaped 
highly resistive sulfide ore body embedded in an average 
continental conductivity configuration excited by a two- 
dimensional, non-uniform source which approximates an 
ionospheric electrojet current. A horizontal dipole above 
a square ore body embedded in poorly conducting rock 
constitutes the third model. Profiles of the electromag- 
netic field quantities showing source effects are 
presented. 

A comparison of theoretical calculations with 
analogue model results is made for uniform two-dimensional 


non-uniform and a horizontal dipole sources over a three- 
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dimensional anomaly. The anomaly on the geophysical scale 
represents an ore body of lower conductivity embedded in 
poorlycéonducting rock.. It is found that the theoretical 
sources used well represent the analogue model sources 
constructed. In spite of the limitations in the two 
methods the results indicate a good degree of compatability 
between the two methods for studying these problems and 
they can be used together to learn more about the effects 
produced by such anomalies. Also, by using two approaches 
Simultaneously, a better understanding of the characteris- 
tics of each approach is obtained. 

From the theoretical calculations the transfer 
functions and related perturbation and induction arrows 
associated with the electromagnetic field solutions are 
computed. The results show that little source effect is 
noted on the induction and perturbation arrows for the 
models chosen. The arrows can be used in the traditional 
manner to indicate the flow of anomalous currents and 
point in the direction of conductivity inhomogeneities. 
The perturbation arrows do not well represent anomalous 
current vectors. However, by combining the p and q arrow 
vectorially the spatial extent of the conductivity anomaly 
can be outlined. The induction and perturbation arrows 
are stable for a variety of anomalous and normal field 


calculations with the perturbation arrow exhibiting the 
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most stability. For small angular spread in the source 
field polarization skewing of the induction arrows is 


produced. 
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CHAPTER | 
INTRODUCTION 


acl General 

External sources which consist of time varying 
Current systems of finite extent flowing in the ionosphere 
and above produce magnetic fields that induce electric cur- 
rents in the earth's interior. These currents in turn 
produce magnetic and electric fields which contribute to the 
total electromagnetic field at the surface of the earth. 

The induced currents are perturbed by variations in the con- 
Guctiyv ity Srerine of the earth and consequently the 
associated electromagnetic fields are modified. The analysis 
of these modifications can yield an estimate of the electrical 
rondgue tivity distribution, of. the earth. 

In this thesis the global and local problems in 
electromagnetic induction are discussed and the three dimen- 
sional local perturbation problem is considered together with 
various methods for its solution. The mathematical formula- 
tions of the electromagnetic field solutions fof uniform, 
general two-dimensional non-uniform and horizontal dipole 
sources over a semi-infinite layered earth with a plane 
boundary are presented. Using these solutions as the boundary 
conditions for an iterative method the solutions over a 


laterally non-uniform earth are obtained. A comparison of 
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theoretically calculated electromagnetic field profiles 
across a three-dimensional anomaly with the profiles of an 
equivalent analogue model for a uniform, a two-dimensional, 
non-uniform and a horizontal dipole source is made. Perturba- 
tion and induction arrows are theoretically computed for a 
variety of source configurations, frequencies and three- 


dimensional anomalous conductivity structures. 


1.2 Global and Local Problems 

The mathematical problems encountered in the study 
of earth currents are of two types. The first type involves 
inducing fields of global dimensions with averaged conducti- 
vities for the earth as a whole; these are known as global 
problems. The second type may or may not involve inducing 
fields of global dimensions but the conductivity values are 


definitely local values; these are known as local problems. 


1.2.1 Global Problems 

In the first group of problems, since we are 
interested in large regions having dimensions comparable to 
those of the earth, spherical polar coordinates are nor- 
mally used. The electrical conductivity o is then treated 
as a smoothed function of the coordinates (r,6,¢) at any 
point within the earth. This function will not take into 
account any local variations in conductivity (which may be 


considerable) but only large scale variations of suitably 
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3 
defined average o for any region. For this problem the 
currents induced in any earth model will depend on the nature 
and distribution of the inducing field. The first step in 
solving such a problem is to express the inducing field in 
terms of spherical harmonics (assumed given). The induction 
effect of each harmonic is then studied separately and the 
result obtained from the summation of all the harmonics. The 
basic theory and method of solution for this type of problem 
is straightforward, though the calculation may be difficult 


and laborious. 


1.2.2 Methods for Solving Global Problems 

One global problem which has been extensively 
Studied is that of the induction, by a varying magnetic field. 
of electric currents in a spherical conductor in which the 
conductivity is a function of r, the distance from the center 
of the earth. Lamb (1883) treated the case for the currents 
induced in a spherical conductor by variation of an external 
magnetic potential. Schuster (1889) applied the method of 
Spherical harmonic analysis, as shown possible by Gauss (1839), 
to the field of the daily magnetic variation and found that 
the major part is of external origin but that a part is also 
produced from within the earth. The induced part of the 
Magnetic field was attributed to electric currents in the 
earth induced by the external magnetic field which he consi- 


dered to be primary. In cooperation with HLamb,Schuster(1889) 
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concluded that the conductivity of the earth was not uniform 
but that the earth consists of layers with the inner layers 
more highly conducting. Chapman (1918) analyzed the average 
field of magnetic storms and showed that the magnetic fieldis 
mainly of external origin with a minor induced part, further 
confirming the work of Schuster (1889). Chapman (1919) 
divided the earth into two concentric layers and gave an 
estimate of the cdnductivity and size of the uniform inner 
core and outer poorly conducting layer. The conductivity 

of the core was determined to be 3.6 x 10°? mho/m and the 
thickness of the surrounding non-conducting layer was 250 km. 
Chapman and Whitehead (1922) showed that the ocean would have 
an appreciable effect on the induced field. Price (1930, 
1931) extended Lamb's solution to aperiodic fieldsand Chapman 
and Price (1930) then used this extension to discuss the 
induced part of the field of magnetic storms. The more 
general case when o is any function of r was studied by 
Lahiri and Price (1939). They developed both formulae and 
methods of calculation for both periodic and aperiodic fields 
when dealing with any distribution of conductivity within the 


sphere r=a that can be expressed in the form: 
Oo = o,(qa/r) for. F<. Ga <a andio (= .0. for, ro > oa ees) 


where ors, ol and r are constants. 
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5 
Lahiri and Price (1939) applied their results using 
the data of Chapman (1919) to investigate the earth's conduc- 
tivity structure to a depth of approximately 800 km. Their 
results showed that the earth's conductivity was uniform and 


> to 10 hho down to a depth of 700 km where 


of a value of 10° 
the conductivity increased to a value of at least 1 mho/m or 
more. 

Using data from the Second Polar Year (1932-1933) 
Hasagawa and Ota (1948) and Benkova (1940) independently de- 
termined that in Chapman's original model the non-conducting 
layer must be 400 km thick and that the conductivity of the 
core must be 5 x frags mho/m. Rikitake (1950), using the same 
data, determined that the earth is composed of a non-conduc- 
ting layer of thickness 400 km, a conducting mantle of 
conductivity 5.0 x 03! mho/m and a core whose conductivity 
must be at least an order of magnitude above that of the 
mantle. 

Due to skin depth considerations, the maximum 
depth of penetration for induced currents is roughly 800- 
100 km, the conductivity distribution of the mantle, which 
extends by seismological evidence from a depth of 700 
km to 2900 km, is undetermined below a depth of 100 km. 
Below this depth the longer period geomagnetic secular 
variation must be used. McDonald (1957) chose a power 


law for the conductivity distribution of the mantle be- 


low 1000 km (conductivities greater than 10 mho/m) and a 
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6 
compendium of Lahiri and Price (1939) conductivity estimates 
for depths less than 800 km (conductivities less than1 mho/m) . 
The intermediate values were smoothed extrapolations from the 
inner and outer portions of the mantle. Yukutake (1959) also 
used the secular variation to estimate the conductivity pro- 
file of the earth. 

Cantwell (1960) gave a more detailed conductivity 
Structure of the upper mantle on the basis of magnetotelluric 
results. A rapid rise in conductivity near a depth of 80 km 


fOr a’ Value,vin excess of 107¢ 


mho/m was inferred by Cantwell] 
(1960). The Cantwell (1960) and McDonald (1957) curves merge 
at 800 km in depth. 

With the introduction of electronic computers, 
time-series analysis of large data sets became possible. 
pe hard’ et al. (1963) studied geomagnetic variations of 
long period and found strong spectral peaks at periods of 
six months and harmonics of 27 days. The data fitted the 
McDonald (1957) profile quite well to depths of approxi- 
mately 1000 km. Banks (1969) extended this method further 
so that not only the spectral peaks but also the continuum 
of the geomagnetic spectrum could be used as well. In 
the model that was obtained a steep rise in the conduc- 
tivity took place at 400 km instead of 700 km as in the 
McDonald (1957) and the Lahiri and Price (1939) models. 
Further analysis by Banks (1972) found that the depth 


of the rise of conductivity is much closer to 700 km. 
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The variation of o with 6 and » within the crust 
is of importance in both global and local problems. For 
example, the conductivity of the oceans is several orders of 
magnitude greater than that of the continents. Since the 
oceans are a thin surface layer on the earth of relatively 
high conductivity, the problem can be modelled in terms of 
electromagnetic induction in thin sheets. Most models consist 
of thin plane or spherical sheets. The conductivity of the 
underlying earth represented by a parallel conducting plane 
Or spherical shell of suitable conductivity at a suitable 
depth. When the effects of the ocean as a whole are con- 
Sidered the plane is infinite or the spherical shell is 
complete. When coastlines are being studied the plate is 
finite and in the spherical case a spherical cap is used. If 
it is desired to model the actual conductivity distribution 
of the real oceans of the earth numerical methods mustbe used. 

The earliest work on induction in thin sheets was 
done by Maxwell (1891) for uniform plane and spherical sheets. 
Lamb (1887a,b) studied non-uniform circular disks. More 
recently Price (1949) developed a general theory of induction 
in thin sheets and shells and obtained boundary conditions at 
the surface of the sheet in terms of the non-uniform distri- 
bution of conductivity and the scalar magnetic potentials of 
the induced and inducing fields. Price also proposed two 


iterative methods for solving the problem. In the first 
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8 
Bailey (1970) proved that if the electromagnetic response of 
the earth is known for all frequencies in a spatial distribu- 
tion represented by a single spherical harmonic the conduc- 
tivity distribution of the earth can be uniquely determined 
assuming spherical symmetry for an exact inversion method. 
Parker (1970) applied the Backus-Gilbert(1970) inversion method to 
Banks' 1969 data and obtained significantly greater surface 
conductivity values than Banks (1969). Parker (1972) obtained 
bounds on the conductivity in the mantle and tested the 
observational data for self consistency. Jady (1974) applied 
an eigenvalue approach to the inversion problem andshowed that 
the top 600 km of mantle is of relatively low conductivity 
with a rise in conductivity taking place somewhere between 
600 and 950 km in depth. Using the data of Malin (1973), 
Jady (1975) determined that low conductivity is required to 
a depth of 650 km below the surface. 

Apparently there have been no successful attempts 
to obtain solutions (analytical or numerical) for the induc- 
tion of currents in a spherical conductor, when the conduc- 
Livi ty.isenor only, aofunction of, r-but also, depends, upon 6 
or ¢ or both. It seems reasonable that, except within the 
outermost 700 km of the earth, the variation of the conduc- 


tivity with r is far more important than any dependence upon 
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9 
method self induction and mutual induction are ignored in 
the starting approximation which is an appropriate approxima- 
tion for a low frequency inducing field. The second method 
appropriate to a high frequency inducing field assumes that 
the vertical component of the induced field is equal and 
Opposite to that-of the inducing field in the starting 
approximation which is equivalent to treating the thin 
conductor as perfectly conducting initially. Bullard and 
Parker (1970) used the first method to investigate the 24- 
hour component of the Sq field in the actual oceans. 
Difficulties are encountered in the second method when 
conductivity discontinuities exist (Hobbs, 1971). However 
a solution in terms of an asymtopic series in inverse powers 
of the frequency can be obtained in some cases (Price and 
Ashour, 1974), Hutson, Kendall and Malin(1972) have reformu- 
lated the first method of Price (1949) so that high frequency 
problems can be solved and proved that for axisymmetric 
problems the solution converges. For non-axisymmetric prob- 
lems no proof exists, however, the method has been success- 
fully applied by Hewson-Browne (1973). Hewson-Browne et al. 
(1973) applied a modified version of the method to the 
problem of induction in a thin strip of finite conductivity 
and obtained a solution which agrees with Parker (1968). 
Parkinson (1975) investigated the possibilities of combining 


the two methods. The first method was used in the ocean 
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10 
regions and the second method was used over the continents. 
Although different iteration algorithms were used the conver- 
gence rate of the problem was still insufficient to give 
useful results. 

Exact solutions can be obtained for only a few 
conductivity distributions which include perfectly conducting 
finite sheets. Smythe (1968) and Lamb (1945) solved the case 
of an infinite strip of uniform width and a circular disk for 
a uniform inducing field using elliptic and oblate spheroidal 
coordinates. For more complex inducing fields, Ashour (1965a) 
solved the problem exactly with a system of dual integral 
equations. The solution for spherical caps subjected to 
symmetric inducing fields has been obtained by Collin (1961) 
and Ashour (1965a,b). Induction in a highly conducting 
hemispherical shell was considered by Doss and Ashour (1971) 
with the purpose of determining the range of frequencies for 
which an ocean of global dimensions acts as a perfect conduc- 
tor except near the coast line. The non-symmetric inducing 
field case can be reduced to a system of integral equations 
(Ashour, 1965a,b). Weidelt (1971) studied a two-dimensional 
electromagnetic induction problem for twoadjacent hal f-sheets 
of different conductivity. The problem was solved by contour 
integration and showed that the current and tangential 
hiddie eroeRnerd component are both finite and discontinuous 


at the boundary of the conductivity discontinuity. A similar 
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result was obtained by Parker (1968) for induction in an 
infinite strip by a parallel line current. For sheets in the 
form of a surface of revolution and when the conductivity 
and inducing field have axial symmetry, the induction problem 
can be solved in terms of a Fredholm integral equation 
(Ashour, 1950). Roden (1964) considered the case of a 
uniformly conducting strip with an infinitely conducting 
plane beneath and reduced the problem to a Fredholm equation. 
A sheet in which the conductivity decreases near the edge to 
zero represents a more realistic model of the oceans. Ashour 
(1971a,b)considered a circular disk, an infinite uniform strip 
and a hemispherical shell with such a conductivity contrast. 

Exact solutions for perfectly conducting closed 
sheets includes the work of Ashour and Price (1948) for 
induction in a spherical shell for a non-uniform conductivity 
and inducing field, Price (1949) for induction im a plane 
Sheet with harmonically varying conductivity and Ashour and 
Ferraro(1964) for applications to the ionosphere for ani- 
Sorropic. conductivity distributions. More recently “Hobbs 
and Price (1970) have derived surface integral formulae for 
expressing any one field quantity in terms of any other for 
currents flowing in concentric spherical surfaces. This 
method is very useful for determining the vertical magnetic 
field component. Hobbs (1971) successfully applied these 
formulae to a hemispheric shell representing the Pac iit ie 


Ocean surrounding a perfectly conducting sphere. 
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For very slowly varying magnetic fields it is not 
PoseVole To treat, the icrus# as a thin. shel J: by. itself, 
because of the mutual induction between the crust and conduc- 
ting. mantle below, ‘Cox.(1960), first_pointed out that the 
conducting mantle should be taken into account when modelling 
electromagnetic effects over the ocean. Rikitake (1961) 
investigated the problem with two-dimensional analog models 
and a theoretical mathematical model of a conductor covered 
with a conducting sheet. Bullard and Parker (1970) proposed 
a model for the conductivity profile of the oceans in which a 
conducting mantle was overlain by a non-conducting crust- 
mantle which in turn was covered with a layer of conducting 
sediments of variable thickness. Cox (1971) assumed that the 
ocean and its sediments were bounded by perfect insulators 


which prevented current leakage into the conducting mantle. 


i-c.o Local Problems 

In the local problem limited regions of the 
earth are concerned and variations of conductivity 
within distances of the order of 100 km are 
studied. For this ‘the sphericity of the earth can be 
ignored and the earth treated as either a semi-infinite or 
thick plate conductor having a non-uniform distribution of 
conductivity. The inducing field in this problem is often 
of global dimensions and therefore effectively uniform over 


the region of interest. However, a knowledge of the inducing 
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field and surrounding ground conductivity in the immediate 
Vicinity of a particular ‘station are insufficient to define 
the strength of the induced currents in that region, as 
pointed out by Price (1965). The induced currents are affec- 
ted by the distribution of the entire inducing field and by 
the average properties of the conductor over a region of 
corresponding dimensions. The entire induced current system 
will contributeto the field at a particular station, but the 
largest contribution comes from induced currents in the 
immediate neighborhood. The important point is that not only 
the local properties, but the properties of the conductor asa 
whole determine the strength of the induced currents in the 
region of interest. 

The unitormity of the tnducing field for the jocal 
problem has sometimes been taken to imply that the actual 
distribution of the inducing field can be ignored when 
calculating the induced currents and fields in the region of 
interest. Ih1S assumption iS impticit in the theory of the 
"magnetotelluric method" of conductivity sounding of Cagniard 
(1953). For models in which thie CONGUCELVILY 1s a funmcr rons oT 


depth only, Cagniard defines the “apparent resistivity" as 


] ak 
= —_———_ ——— Liz2 


in which w is the angular frequency of the horizontal, 


orthogonal electromagnetic field components EY and Hy and Uo 
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is the magnetic permeability of free space. The model is 
assumed to be horizontally stratified in the x- and y- 
directions extending to infinity. The formula is derived 
from the skin depth effect for uniform half space conductors 
excited by alternating currents. As such no source field 
effects are taken into account. As w decreases the currents 
penetrate deeper into the conductor and the manner in which 
P, varies gives an indication of the conductivity profile. 
The accuracy of the profile is dependent on the particular 
geophysical situation being investigated. For moderately 
uniform terrain and depths the method is convenient and valu- 
able. At greater depth even for horizontally stratified 
strata it may be necessary to include consideration of the 
source field as discussed by Wait (1954) and Price (1962). 

Price (1964) points out that the variation fields 
Should first be separated into parts of external and internal 
origin and examined separately in order to objectively 
interpret local patterns of geomagnetic variations. The 
anomalous patterns in the fields of the induced currents are 
produced by non-uniform distributions of conductivity in the 
earth, since the external inducing field can be assumed to be 
of a relatively simple form and varies only slightly over the 
anomalous area except near the equatorial and auroral elec- 
trojets. There seems to be no reason to suggest that the 


ionosphere has anomalous properties over a specific geogra- 
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15 
phic region, even though a small effect on the (anisotropic) 
conductivity of the ionosphere might be produced by a strong 
perturbation, of .the earth's, magnetic field....Even if a 
considerable disturbance is caused in the ionosphere by 
such local anomalies it is doubtful that the surface values 
of the magnetic field would be greatly affected. Also, if 
the external field is found to reflect to a small extent the 
anomalous magnetic field distribution, this can be attributed 
to the general mutual induction between the earth and the 
ionosphere. 

Once the external field and its induced current 
system is known, the local problem can be reduced to examin- 
ing the local redistribution of a given average system of 
induced currents caused by local inhomogeneities of conduc- 
tivity. However, it is not necessarily a simple problem of 
redistribution of steady current flow, since the currents are 
time varying and therefore a skin effect will be encountered. 
The mathematical problems which must be considered are those 
relating to the disturbance of skin effect distributions ofa 
current system rather than actual electromagnetic induction 
problems in which the distribution of the induced currents 
for a non-uniform conductor is determined directly from a 
given inducing field. This problem is described by Price 


(1964) as follows: 
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"Using cartesian coordinates (x,y,z) 

with the z axis vertically downwards, a non- 

uniform conductor occupies the half-space 

74720. aNear, thesorigin the conductivity. dis 

a function (not necessarily continuous) of 

(xXoyz). but at great. distances, from. the 

OVAGin dts is,.as function. ofazeoniy .us Acgiven 

alternating e.m.f. impels currents near the 

Surface of the conductor. The problem is 

to determine the distribution and surface 

field of these currents." 
The given e.m.f. is the electric field arising from the 
varying magnetic field and is controlled by the properties 
of the conductor at great distances. In essence the problem 
is formulated in terms of the redistribution of a known 
current system near a non-uniform conductor in which the 


distribution of the current system is known at large x and y. 


1.2.4 Methods for Solving the Local Induction Problem 

In his classic paper, Price (1950) considered 
electromagnetic induction in a semi-infinite conductor with 
a plane boundary for any inducing field, which was assumed 
known. A uniformly layered earth model was considered by 
Tikhonov (1950). . Cagniard (1953), in his derivation of the 
well known theory of magnetotellurics, assumed a spatially 
uniform source of infinite extent and sinusoidal time 
variation, as well as a semi-infinite uniformly stratified 
earth with a plane boundary. Wait (1954) and Price (1962) 
discussed Cagniard's results in terms of a source of finite 
dimensions. However, until the last decade, little attention 


has been paid to the problem of a vertical fault in the 
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conductive region which represents the earth. It is well 
known that there are limited regions in various parts of the 
earth where temporal changes in the geomagnetic field at 
relatively near stations show considerable differences in 
form and amplitude. These differences remain consistent in 
Character for the fluctuations of a given frequency, and 
form a definite pattern for the area studied. They are 
associated with lateral variations in conductivity below the 
earth's surface. Two-dimensional magnetic variometer array 
studies have well exhibited the consistency of the field 
configurations and that the differences in form and amplitude 
are associated with conductivity variations. It is therefore 
important to understand the manner in which lateral varia- 
tions affect the surface field in order to interpret such 
variations. 

D'Erceville and Kunetz (1962) obtained a solution 
for two media of different conductivities in contact along a 
vertical plane overlaying a horizontal basement that iseither 
infinitely conducting, infinitely resistive-or at infinite 
depth when the magnetic field is everywhere parallel to the 
strike of the fault, (H-polarization, see Chapter 2). Rankin 
(1962) applied the method of d'Erceville and Kunetz to the 
case of a dike of infinite length but finite depth and in 


which the magnetic field is parallel to the dike. 
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Weaver (1963) considered a half-space conductor 
which consisted of two quarter-spaces of different conduc- 
tivity with a plane vertical fault of infinite depth. He 
considered both the H-polarization case and the E-polariza- 
tion case (see Chapter 2) using an analytical technique. 

His H-polarization case agreed well with that of D'Erceville 
and Kunetz. However, in order for him to obtain a solution 
for the E-polarization eee he found it necessary to use the 
approximate boundary condition that the tangential magnetic 
field is constant along the surface of the conductor. Mann 
(1970) proposed a perturbation technique to consider the 
original Weaver (1963) problem, and has shown that 

Weaver's original solution is just a first approximation to 
the field inside the conductor for normally incident 

waves. Weaver and Thomson (1972) have applied this tech- 
nique and have been able to avoid the approximate boundary 
condition of the constancy of the horizontal magnetic 

field along the surface and have obtained approximate 
solutions. An approximate solution to the field was found 
for a periodic line current above a non-uniform earth 

when the height and magnitude of the line current approached 
infinity in such a way that the inducing field near the 
earth became uniform and finite. Geyer (1970) has used a 
similar perturbation technique to investigate the electromag- 
netic anomalies over several types of subsurface structures, 


and more recently (Geyer,1972) has extended this to a dipping 
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contact. Previous to this, Yukutake (1967) considered 
induction in a conductor bounded by an inclined interface 
with a small angle of tilt by a successive approximation 
method which used repeated reflections of electromagnetic 
energy between the ground surface and the tilted boundary. 
Treumann (1970a,b,c) considered induction in non-uniform 
plates of finite thickness, and has been able to obtain 
approximate solutions for the field at the surface of the . 
plates when the external inducing field is uniform. Weidelt 
(1971) has studied induction in two adjacent half sheets with 
different uniform conductivities. Furthermore, Schumucker 
(1971a) has used convolution integrals to investigate the 
induction in a model with a non-uniform surface layer above a 
layered substratum. Hughes (1973, 1974) presented a semi- 
analytical linear solution to the two-dimensional E-polariza- 
tion case problem involving Fourier series representations of 
the conductivity anomaly. A two-layered earth with a 
Sinusoidally varying overburden was modelled by Hughes and 
Wait (1975) using the semi-analytical approach of Mann(1964). 

Analytical methods generally give solutions for 
only specialized cases. Anomalous conductivity structures 
of any shape may occur naturally, and so it is necessary to 
consider methods to deal with inhomogeneities of arbitrary 
shape. Several numerical methods have been used to obtain 


solutions for problems with vertical discontinuities for both 
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H-polarization and E-polarization cases. Dulaney and Madden 
(1962) were the first to apply the transmission line analogy 
over a two-dimensional mesh to such problems as these. The 
transmission line analogy arises from the similarity in form 
between Maxwell's equations governing the orthogonal compo- 
nents of E and H and the transmission line equations governing 
current and voltage on a transmission line. This method has 
been used by many authors (including Madden and Thompson, 
meo, Madeaen and*Swift, 1969; Swift, 1967, 1971; Wright; 
mote 0/0, slankis., 19705) VozottT, 19/1 ey Kku-et ats. 1973. 
and Ku, 1976). The finite element method, which uses the 
principle that electromagnetic fields behave in such a way as 
to minimize the energy, has recently been applied toinduction 
problems by Coggon (1971), Ryu (1972), Silvestor and Haslam 
(1972), and Reddy and Rankin (1973). An anisotropically 
conducting laterally inhomogeneous mediumwas modelled by 
Reddy and Rankin (1975) using the finite element method. 
Rodi (1976) developed a finite element method which increases 
the computational accuracy of the field components used in 
hagnetote tLuric 4 calculations.. “Finite difrerence cechniques 
have been applied by Neves (1957), Latka (1966), Patrick and 
Bostick (1969) and others. Jones and Price’ (1969, ° 1970), 
1971a,b) have employed the finite difference technique for 
studies of various two-dimensional conductivity distributions, 


and Jones and Pascoe (1971) and Pascoe and Jones (1972) have 
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21 
given a general computer program for the solution of the 
local perturbation problem. 

The relative merits and disadvantages of the four 
numerical techniques used to solve the local electromagnetic 
induction problem can be used to decide which method best 
f7tS a parbicular, mode) ling, Situation... Both, the. finite 
element and transmission line analogy methods use matrix 
inversion to obtain solutions. For properly chosen grid 
Spacings economical solutions in terms of computer costs can 
be obtained. However, grids with large numbers of points 
cannot be accommodated as the computer storage for the 
matrix inversion becomes prohibitive. Finite difference 
techniques when properly applied can be used with large 
grids since matrices are not involved. The computing time 
used is somewhat longer and the stability of the solution 
tends to be less than for the finite element method solution. 
The integral equation method reduces the domain of interest 
to the anomalous conductivity structure, and the electro- 
magnetic field solutions can be obtained in terms of the 
scattering currents within it. However large anomalous 
domains cannot be handled as matrix inversion is involved 
and large conductivity contrasts (> 100:1) cannot be 


accommodated. 


The grid mesh size in all the methods must be 


chosen properly to obtain a rapid solution. The smaller 


and more even the mesh size is across conductivity contrasts 
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the more rapid is the convergence of the solution. If 
coarser grids are chosen then the convergence of the 
solution takes substantially longer. However, the total 
number of grid points is limited and therefore the mesh 
sizes are chosen so that small uniform steps are taken 
across conductivity contrasts. Progressively larger grid 
spacings are used as one approaches the boundary of the 
mesh so as to keep the anomaly as far away from the 
boundary as possible. 

Jones and Price (1969, 1970) considered the half- 
Space conductor as the limit of a spherical conductor as the 
radius becomes infinite. They then carefully developed the 
boundary conditions and were able to obtain solutions for both 
the H-polarization and E-polarization cases. Jones and Price 
(1970) studied the field distribution within the whole two- 
dimensional region of interest, as well as the surface values 
of the various components, including the Cagniard (1953) 
apparent resistivity. In more recent work, Jones and Price 
(1971a,b) and Jones (1971a,b) as well as Hibbs and Jones 
(1972 ) have calculated apparent resistivity values along the 
surface of the conductor. Tatrallyay and Jones (1974a,b) 
applied a finite difference technique to a cylindrical mesh. 
Losecke and Muller (1975) investigated a two-dimensional 
overhanging, high resistivity structure utilizing a finite 
difference technique. An increase in the rate of convergence 


of the above problem was obtained by Muller and Losecke (1975). 
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Williamson, Hewlett and Tammemagi (1974) discovered 
an error in the Jones and Pascoe (1971) representation of 
the first derivatives for variable grid spacings. Jones and 
Thomson (1974) showed that its effect on the calculated 
results is much reduced when the grid spacings used are not 
too irregular. Brewitt-Taylor and Weaver (1976) show that 
for irregular grid spacings across regions of conductivity 
change the equations for the E-polarization case must be 
further modified to take into account weighted conductivity 
averages. The B-polarization equations were found to be 
incorrect and must be reformulated for all but regions of 
uniform conductivity. Brewitt-Taylor and Weaver regard 
conductivity changes as roughly linearly from one value to 
the next thus reinforcing the Lines and Jones (1973) 
concept of transition zones. 

Tt should be noted that in all the’ numerical 
methods mentioned above the assumption has been made eae 
the inducing field is uniform. A non-uniform source in the 
form of a line current has been considered by Wait (1962) 
and Dosso and Jacobs (1968). Schmucker (1971b) has considered 


non-uniform sources over a laterally inhomogeneous earth in 
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which the lateral changes in conductivity are confined to a 
limited depth range. Hermance and Peltier (1970) studied the 
magnetotelluric field of a line current over a layered 
Subsurface, and later (Peltier and Hermance, 1971) extended 
their approach to include symmetric current distributions, 
such as a Gaussian current intensity source. Hutton (1969, 1971) 
and Hutton and Leggeat (1971) investigated the problem of 
electromagnetic induction by an equatorial electrojet taking 
into account the range of periods over which the electrojet 
effect is exhibited and the relative contributions of the 
"normal" and "enhanced" part of the inducing field. Hohmann 
(1971) investigated the two-dimensional line current source 
problem by deriving a Green's scalar function from Maxwell's 
equations and expressing the electric field in terms of an 
integral equation involving volume integrations of the 
anomalous field vector over only the anomalous domain. The 
boundary conditions were incorporated in the kernel of the 
integral equation and were thus automatically satisfied by 
the solution. The two-dimensional electromagnetic response 
of a buried cylindrical inhomogeneity excited by a line 
current source was calculated by Howard (1972) using a 
similar method. Schmucker (1971b) has considered models 
with a non-uniform source over a laterally inhomogeneous 
earth in which the lateral changes in conductivity are 


confined to a limited depth range. Hibbs and Jones (1973a) 
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have used the method of Peltier and Hermance (19/7 laeto 


obtain boundary values for the method of Jones and Price 
(1970) and Jones and Pascoe (1971) so that perturbations 
of such fields by embedded inhomogeneities in the E- 
polarization case may be computed. Hibbs and Jones (1973b) 
extended this method to consider a non-symmetric, non- 
uniform source, as well as aperiodic spatially time-varying 
sources (Hibbs and Jones, 1974). Rikitake (1975) investi- 
gated electromagnetic induction in the earth for a modula- 
ted inducing field and showed the possibility of appreciable 
error in electrical conductivity estimation due to the 
contamination of long period variations in the geomagnetic 
field by short period variations. The electromagnetic 
response of a two-dimensional, non-uniform earth to an 
oscillating magnetic dipole source was solved numerically 
by Stoyer and Greenfield (1976) using a transmission line 
analogy. Hibbs and Jones (1976a,b) have given a method 
for modelling any arbitrary two-dimensional source field. 
The use of two-dimensional arrays of magnetic field 
measuring instruments covering large areas of the earth 
such as the Western United States and Canada (Porath et al., 
io ies Reat ze) etical.:, W9/0s Canpheld. et at ao. ‘Pore Gh 
et: ates 197 le Porath «and Gough, gUo7 le -Camiield. + lara. 
Camfield and Gough, 1975; and Alabi et. at... 1975), South- 
eastern Australia (Bennett and Lilley, 1972; Lilley and 
Bennett. 1972: Gough: et ah... 9/9 Bennett and Lilley, 


1973; Bennett and Lilley, 1974; and Lilley, 1976) and 
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Africa (Gough et ali, 1973:°and.de'Beer et alsiive 9.97 6) 


and studies related to the interpretation of array data 
(Lilley and Tammemagi, 1972; Tammemagi and Lilley, 1973; 
Gough, 1973; Frazer, 1974; Garland, 1975; and Lilley,1975) 
have stimulated interest in the three-dimensional induction 
problem. Treumann (1970d) first showed that the Green's 
tensor could be used to obtain a solution to the three- 
dimensional problem. Jones and Pascoe (1972), in an exten- 
sion of their two-dimensional finite difference method, 
obtained a solution for a three-dimensional buried conduc- 
tivity inhomogeneity using a cubic mesh. This work was 
extended to grids of variable dimensions by Lines and Jones 
(1973). Hohmann (1975) extended his two-dimensional model 
theoretical solution to three-dimensions involving a dyadic 
Green's function. The theoretical solution was obtained 

in terms of an integral equation which was reduced to a 
matrix equation numerically soluble for the scattering 
current in the anomalous body. The electric and magnetic 
fields outside the body were found by integrating the 
appropriate half-space dyadic Green's function over the 
volume of the scattering current. Raiche (1974) and Raiche 
and Coggon (1975) approached the three-dimensional problem 
in the same manner. Howard (1975) studied the problem of 

a buried cylinder excited by a magnetic dipole in three- 
dimensions. Also, Weidelt (1975) used the Green's tensor 


method to calculate the electromagnetic induction effects 
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of a three-dimensional structure. Hibbs and Jones (1976c) 
have calculated the electromagnetic field response for 
various three-dimensional conductivity structures for a 
uniform, two-dimensional non-uniform and a horizontal 


dipole source. 


12-2 “Transfer Functions 

The: use of transfer functions for the:.detection 
of anomalous internal currents originated in the work of 
Parkinson (1959) and was further advanced by Wiese (1962). 
Parkinson (1959) observed an empirical relationship between 


the temporal variation of the horizontal (AX, ,AY and 


1) 
vertical (AZ_) magnetic field components. This linear 
relationship 


SOA AK ot BAY 


r t (ele) 


AZy 


where A and B are constants to be determined, has come to 
be known as the "Parkinson Relationship". Schmucker (1970) 
extended the original work and provided a formal theoretical 
approach to the description of the anomalous fields in 
terms of transfer functions and described both induction 
(or Parkinson arrows) as well as perturbation arrows. 
Everett and Hyndman (1967) introduced a unit vector method 
of transfer function estimation and Cochrane and Hyndman 
(1970) and Hyndman and Cochrane (1971) simplified the cal- 
culation by substituting the observed field components in 


place of the normal field components in the calculations. 
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Edwards et al. (1971) presented another approach to the 
transfer function calculation in which they considered 
magnetic correlation (M) vectors, and plotted Parkinson 
erreows forystations in the®British@istessvrhHonkura ¢1971) 
compared experimentally computed induction arrows using 
data from Miyake-jima island with a theoretical model 
calculation of the arrows. Lilley (1974) analyzed the 
mathematical formalism of Schmucker (1970) in terms of an 
“induction tensor", and compared its characteristics for 
two- and three-dimensional cases. The time variation of 
induction arrows was studied by Anderson (1975) for the 
ULC... frequency band (1-100 mHz). « Serious distortions in 
the inferred induction arrow directions due to source field 
variations were shown by Anderson et al. (1976a,b). 

Weidelt (1975) calculated induction arrows for 
a theoretical model of an embedded anomaly using two 
mutually perpendicular polarizations of the external elec- 
tromagnetic field to construct his arrows. Hibbs and 
Jones (1976d,e) have calculated theoretical induction and 
perturbation arrows for uniform and two-dimensional non- 
uniform source fields oriented in various directions 
relative to an embedded three-dimensional anomaly. A 
similar calculation has been made (Hibbs and Jones, 1976f) 


for a horizontal dipole source over a highly conducting 


anomaly for two dipole positions. 
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The work of this thesis is primarily concerned 


with finding solutions to the three-dimensional local 
induction problem for various source configurations and 
determining the source effect for each source type on 
the transfer functions which are calculated from these 
solutions. A comparison of theoretical solutions with 
an analogue model equivalent for different source field 


types is-also presented. 
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CHAPTER: 2 


SOLUTION OF THE THREE-DIMENSIONAL LOCAL INDUCTION 
PROBLEM FOR ONE, TWO AND THREE DIMENSIONAL SOURCE FIELDS 


In this chapter a method of solution for the three- 
dimensional local electromagnetic induction problem for 
laterally non-uniform conductivity distributions is given. 
The general local induction problem as described by Price 
(1950) is discussed. The solutions for three source config- 
urations (uniform, two dimensional non-uniform and a horizon- 
tal dipole) above a uniform layered conductor re presented. 
In the method, these solutions are used to give boundary 
conditions for the Lines and Jones (1973) finite difference 
technique from which the solutions for the electromagnetic 
fields associated with three-dimensional conductivity 
inhomogeneities are obtained. The conductivity structures, 
which differ for each source configuration used, are chosen 
to represent varying degrees of symmetry and conductivity 
contrasts. For each combination the source frequency is 
chosen such that the amplitudes of the electromagnetic field 


components presented will have large values. 


2.1 The General Three-Dimensional Local Electromagnetic 


Induction Problem 


Price (1950) treated the general theory of the 


induction and of free decay of electric currents in a semi- 
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infinite conductor with a plane boundary. Two fundamental 
types of solutions were found and physical explanations for 
these were presented by him. Methods of calculating induced 
fields and their current systems corresponding to external 
inducing magnetic fields were obtained. 

The mathematical development here follows closely 
that given by Price (1950), with the exception that the M.K.S. 
system of units is used throughout. For a continuous medium 
the behavior of the electromagnetic fields is described by 


Maxwell's equations 


aD 
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Together with the relations 
ee Gee se I er D = p : (2.<3') 


where 0 is the volume distribution of charge. In the 


isotropic and uniform media which are considered throughout 


this work 
BS eh ae 0. Fe Oe Oe tee (2.4) 


where the permeability, uy, the conductivity, o, and the 


permittivity, ¢, are constants. 


Operating on both sides of equation (2.1) with the 


divergence operator and using the appropriate substitution 
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and by using (2.3) we have 
o op _ 
SC) cae 0 (2-263) 
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0 
-=t 
Boe ces ; ener) 


which shows that any initial volume charge distribution (p,) 
in the conductor disperses at a rate which is not influenced 
by other field conditions. Therefore, inside a uniform 
conductor the volume charge can be neglected without loss of 
generality if the relaxation time constant satisfies the 


condition 
<<] 226) 


We have, then, 
Ve Foca an! ( 2u Qa 


By taking the curl of (2.1) and substituting the expressions 
Oorained for ¥ x E into (2.2) or piermecivetye oe taking the 
curl of (2.2) and substituting the expression obtained 

for V x H into (2.1), it can be shown that the electro- 
magnetic field quantities satisfy an equation of the 


form 
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2 
2 e) 0 
V G = U 8 aby BS (2 .10.) 
t 
where G’may*be either E or H. In a conductor, when the time 


changes in the field are not too rapid, we 3G/at can be 
ignored when compared to woG. iat is, displacement currents 
are negligible when compared to the conduction currents in 
mesconductor. Garland (1971) discusses the conditions for 
which displacement currents may be neglected. From equation 
(2.1), if L and T are taken as the characteristic length and 
time of the variation of the electromagnetic field, aD/ot is 
bry tnerorderse E/P andi VcxrEciswofhtnhnetordergkf~h soaspdl/te 
I< wollowsothatmoD/ dtotsaofethetorder peml Eater Also V x H 
is of the order H/L. For oD/at to be negligible in compari- 


son with V x H 
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Hence, the characteristic time must be great compared to the 
time of travel of the electromagnetic wave across the region. 
For a characteristic length equal to the diameter of the 
earth the characteristic time must be much less than 0.03 
seconds. In all problems considered in this work the above 


the condition is satisfied. Equation (2.10) then reduces to 


the diffusion equation 
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in the conducting region. 
In a dielectric o is negligibly small and equation 


(2.10) becomes 


C2) 


where the electromagnetic effects propagate with speed 


ylfe 


(ue through the region. Induction fields vary in a 
relatively slow manner, so that the time of travel across 
any portion of the dielectric, within the region considered, 
is negligibly small compared with the time taken for field 


changes to become effective at the point of interest. Hence, 


for a dielectric, equation (2.14) reduces to 
Von G =u 0 (2415) 


The above approximations are equivalent to ignoring 
the magnetic effects of the displacement current in both the 
conductor and the dielectric. Although tne magnetic effects 
of 3D/dt are negligible, this does not imply that the physi- 
cal effects of the displacement current can be ignored. 


Since a ane! is zero, H can be approximately derived from a 


scalar potential, 
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Where & satisfies Laplace's equation. Equation (2.3) shows 
that B is everywhere solenoidal and hence derivable from a 


vector potential A defined by 
Boo ae a A (aah h} 


Substituting this into (2.2) and using the vector identity 


wx (-Vo) ='0, gives 
E=--+3-Vo. . (2218) 


A and ¢ satisfy the fundamental equation (2.10) if they 


satisfy the relation. 
Ty) Ae oie Oe at usd = 0 Carat 9) 


Although the electromagnetic field vectors E and 
B can be determined when A and 6 are known, the functions A 
and o are not uniquely defined since the electromagnetic 
fields are invariant under the gauge transformation. However, 
in a particular problem, any solution of A and ® can be used 
to find the unique fields which satisfy the prescribed 
boundary conditions. 

The electromagnetic field quantities must also 
satisfy certain boundary conditions. The tangential compo- 
nents of E and H and the normal component of B are continuous 
across the interface between any two media. The normal 
components of D and J are not in general continuous and the 


sum of their components normal to the interface are p and 


- ee » respectively. 
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2.2 Elementary Solutions of Maxwell's Equations for a 
Semi-Infinite Conductor With a Plane Boundary 

ihe coordinate system is shown in Fig. 2.1. ° The 
region -|h| < z< 0 is a free space region with conductivity 
zero. A conductor occupies the half-space z > O and it is 
assumed that a current source of some form is located at or 
above z = -|h| 

Since the tangential components of E' and Hare 
continuous at the surface, z=0, and the conductivity, o, 
and permeability, u, are not, the elementary solutions of E 
and H can be expressed in the form Z(z,t) F(x,y) . The 
variable t occurs only in Z because ee is multiplied by uo 


at 
in equation (2.13). Substituting 


sae AG Te) mee ea ®, | (2.20) 


mato. (2.9), (2.10)and (2.15) 
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Since the coefficientsof F on the right hand of 
(2.22) and (2.23) areindependent of x or y they can be 
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represented by a real constant, -A’, so that 
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Figure 2.1 The coordinate system. xX-Y plane represents the 
surface of the semi-infinite half-space conduc- 
tor with a plane boundary. 
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aig ) ue 
ee ee ee i (2.24) 
Xx y 


It will be seen from equation (2.23) that the constant is 
necessarily negative. 


For the conductor 


ee es sé = 0 (2.25) 
dZ 
and for the dielectric 
2 
en (2.26) 
0Z 
Erom (2.21). it follows that’ either 
OF OF 
Gaeta 0 and Crea we = 0 Ade 
oF OF 
ea eer ee 
aH (a5 + a ) ai ae + a (2.28) 


where a is any real or complex constant. 

The elementary solutions corresponding to (2.27) 
and (2.28) are of two types: tne first type corresponds to 
Current systems in the conductor which have magnetic fields 
outside and the second type has magnetic fields which are 


entirely confined within the conductor. 


Ge2.i velementary Solutions of the First Type 


From (2.24) it is seen that F may be written in 


the form 
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F=(22,-2) ,0) | (2.29) 


where P is a function of x and y which from (2.24) must 


Satisfy the relation 


aR < Op 
oe + oF + ep = 0 (2.30) 
0X oy 


From (2.26) the form of Z in the dielectric must be 


2.2 0h then SusuB (Clem e 
so that 
mee {Al te! 4 Be ies , - 2,0) for z<0 = (2.31) 


Since the tangential components of E are continuous at the 


Surface, z=0, E inside the conductor is of the form 
E = 2(z,t)( , - 22,0) for z>0 (232) 
where ZEotOot) = th tt) et abt b) 5 0233.3) 


The magnetic fields inside and outside the conduc- 


tor can be found by using equation (2.2) which gives 
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The field intensity in the conductor must tend to zero as 


z becomes large, so that 
Um I) or Weg Gee ge pe anes C274) 


Ime the region -|h| < z < 0 the potential of the total magne- 
tic field satisfies Laplace's equation and the elementary 


Solution is of the form 


, B'ltte -) Poa) (2.38) 


where A is real and positive (or zero) and P satisfies 
equation (2.30). 

The part of & involving Bee corresponds to the 
source field over the region z < -|h|, that is the 


AZ at 
corresponds to tne 


inducing field. The part involving a 
field arising from within the conductor, that is the inacuced 
field. Therefore A'(t) in (2.38) can be regarded as a known 


function of the source. Comparing with equation (2.34) gives 


ee ee a ee ce (2.39) 
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and A(t) can be obtained. 

EGuations! (233) 90(20 3690 and! C2. 37) together with 
(2625) determine B(t) and Z(zst)}suniquely in terms Of Ar(7t). 
The induced field will be completely determinedif the ini- 
tial conditions, B(-0), are known. The induced field is thus 


expressed entirely in terms of solutions of the first type. 
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In these solutions the vertical component of E and therefore 
the vertical current is everywhere zero. The induced 


currents flow parallel to the surface of the conductor. 


eeco-c | Elementary Soluttons of the Second Type 
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In-this case, inside the conductor from (2.25) 
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where P, is a function of x and y which satisfies (2.30). 
The terms involving Py lead to ‘solutions of the? Hirst type 
and can be omitted. Fy also satisfies (2.30). Hence, by 


letting F_ = AP a second type of solution for E of the form 
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is found. 


Outside the conductor Z satisfies (2.26) and 


comparing with (2.40) it is apparent that a = + A for z< 0. 


* Ta 
snotonon bap 3. to jianoqmo>, pao tetav a sd 028 
‘beoubni dT .ovss svoriwy rave ae) . ri she 


-YOtouUbHOD SHY To. 996TH SAT os otto 


aah brosee Sait. tO. Bur ch (i etname — 


to od Soum XS noiy bettetdse er (as. $) ci WS | 


(O8.S) | shsiyyo = 
| f : 
(28.8) most yotoubnoo 3% sbiznt ,seno BF 


(F®.S) (ap6 big = PVE = 8) 
($8.8) gy 5) tap 
(OK.S) x ala a Bini 
. My 
De | 


sant cae Sd>m69 th (RS. s); wae” 
ace. ML ee | pit, OME ara, 


fee aoe ee ET ter) 2 eee aioe 
(AR.S) i e g, ‘i, ao } eR 


4 
- 


vi 
“tO. 9) esttetise fotty y aie Ki ate: 2t 59 

2 ae 

sqyt tevit Sad Fo enotsuloe’ of bel Ae sia 

yd .sonsH .(0€-S) so Hetdee o2te. ie 


(28.8) Dees ioe ee cue aie 
iwiAs j <1 te 


[ae “bn, (a8. 8) eattaitse 3 v0d3u 


. +S iia ie 
o> Ss vot i + = eas ai a r (OF 


ee ae (Th, a eee 


= »), ue ca - 2 Te | - “ay me 
ee a 


42 


Egis7then,ofs the. form 


TAZ 


E= V{A,(t)e “7 + By(t)e’*2} P(x,y) for z<0 . (2.46) 


The corresponding magnetic field inside the conductor is of 


the form 
OBr i « sBt=az Z oP. dP 
vi Cie (A - @ DV a Dh giaeey O°)F Fa Otome Zee as) 


Sincey ke ise-avec tor outside, thejiconductor 
cn tan OOP Sete yerec0 : (2.48) 


Therefore, solutions of the second type have no 


magnetic field outside the conductor and solutions of this 


type do not enter into the problem of finding currents in- 
duced by external magnetic fields. Solutions of this type 
correspond to the decay of current systems of a- special 
form in the conductor. Freely decaying current systems 
corresponding to solutions of the first type can also be 
obtained by Setting A(t) to zero in (2.38). The modes of 
free decay are discussed in detail by Price (1950). 

By using solutions of the Tirst type, tne? loca! 
induction problem may be solved for uniform, two-dimensional, 
non-uniform and horizontal dipole sources over laterally 
homogeneous conductors. Solutions for the first two of 


these sources can be formulated in terms of the two-dimen- 


sional induction problem. The horizontal dipole source must 


be dealt with in a more general manner. 
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2.3 The Two-Dimensional Problem; E and H Polarizations 
The coordinate system is as before (Fig. 2.1) with 
the region -|h| < z < 0 a free space region withconductivity 
equal to zero, and the conductor occupies the half-space 
Zee eo. Again, 1t 1S. assumed that a current source of some 
form 1s located at‘or above z = =|h| 
For an oscillating field with time dependence 
aaek where w is the angular frequency and the period 21/w 


is sufficiently long so that displacement currents may be 


neglected, Maxwell's equations become: 


curl “i (2.49) 
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where o is the conductivity appropriate for each region, 
Uo is the magnetic permeability of free space and the time 
factor q (iwt) is understood in all field quantities. 

In the two-dimensional problem the fields do not 
wapy in one direction, in this case the x-direction. That 1s 


all quantities are independent of x, and Maxwell's equations 


reduce to 
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JES JE 
yar 55 = Twu Hy 3 (2.54) 
JF 
he Meets (2.55) 
JE ; 

2 ays Le TWH OHS . (2756) 


These equations are such that only EL MN, and He 
Poem vo lived win (2.51..2.55 and 2.56.) and ony Hy EY and ES 
are-involved in (2.52, .2.53 and 2.54)... These two separate 
sets of equations can therefore be solved independently. 


fombi ning, (2.51) .. (2.55) and. (2.56). we obtain 


aE, lee 
—y~ + —5* = iw oE (2557) 
ay a2 ge 


which is the equation to be solved in all regions for the 
E-polarization case (the case in which the eysetene tTretd 1s 
parallel to any conductivity discontinuities). 

Combining (2.52) (2.53) and (2054) a simitar 
equation is obtained for the H-polarization case (the case 
in which the magnetic field is parallel to any conductivity 


discontinuities): 


Ae, 5 = Tw oH, ; (2.58) 
dy ny: 


With the proper values of co inserted for each 
conductive region, and with the proper boundary conditions, 
these equations can be solved for the two separate cases to 


obtain the field distributions. For the E-polarization, EY. 
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is determined from the solution of (2.57), and then (2.55) 
and (2.56) may be used to calculate the two non-zero magnetic 
field components. associated with this case. Also, for. the 
H-polarization, hin may be determined from (2.58), and 
equations (2.52) and (2.53) may be used to obtain the 


appropriate electric field components. 


2.3.1 Determinacy of the Two-Dimensional Problem 


Not all two-dimensional problems are determinate. 
It is important to examine the conditions for determinacy 


for a periodic inducing field over a semi-infinite uniform 
conductor with a plane boundary occupying the half-space 


z > 0. Following Price (1950), we let ES be of the form 
E = Zi zieak Wye ew (2.59) 


and substitution of this expression for EY into equation 


(2.57). gives: 
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for z < 0, where the arguments for Y and Z have now been 


dropped. 
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Since the coefficients of Y on the right hand side 
of (2.60) and (2.57) are independent of x and y they must be 
constants, eae where A is real, positive and non-zero such 


Ghat 


2 
Sed i (2.62) 


dy 


Therefore, in the region z > 0 


2 


Oe Mo 4 8 Za (2,625 
QZ ° 
and ee er Re (2.64) 


where 92 = TWU Yo + 2 and A and B are constants. In the 


regron z.< 0, 


See! i Ye he, (2.65) 


and Z(z) = ae + be : (2.66) 


where a and b are constants. The corresponding E field is 


then 
-0Z 4 
EY = {Ae + Be HY (2 abu) 
fee O wand (oF. 5 fa Paes ae (2.68) 
foe 2 a 0 


As z becomes large, EY must tend to zero so that 


B must be zero. At the boundary z=0, EY must be continuous, 


ae eat a : (2.69) 
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The corresponding magnetic fields are obtained from (2.55) 


mde. 56). for z-> 0: 


Hy = ge gph - Whe ey (2.70) 
Wt gh ges eee Pe eae 
ane for z <, 0: 
Hy = - i Pe ee a en § Cees) 
Me = - ae Gat eae heae esr ony (2th 3) 


The tangential component of the magnetic field must 
be continuous across the boundary z=0, which leads to the 


relation 


eA = r(a-b) . (2.74) 


The scalar potential 2 of the total magnetic field 
in the region -|h| < z < 0 satisfies Laplace's equation and 
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6 = bit ee oo ee mg) (2.75) 


‘ . 2 a 
where P(y) satisfies the equation 9 P(y)/oy = - r*P(y) 


Since the conductivity outside the conductor is 


Ze: O 
vk He © (2.76) 


and HY and o can be found from the relation 
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Pomoanming (2.69), (2.74) and (2.80) gives: 


p= (ei+c , (2.81) 
a = vag (24) 0, (2.82) 
b= - (24a , (2.83) 
and a= - oo, (2.84) 


Therefore by specifying C all the other quantities 
will be determined. This analysis, while confined to uni- 
form subsurfaces, can be extended to multi-layered cases. 

The part of 2 in (2.75) which involves anu 
corresponds to the inducing field of the current sources and 
the part which involves ar corresponds to the field of the 
induced currents. It can be assumed that the inducing field 
potential is specified by the current distribution so that 
C a AZ - P(y) is known. From the preceding discussion it can 


‘be seen that all other fields and potentials are determined 
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from the inducing field. Therefore the problem is deter- 
minate for A positive and non-zero. The H-polarization case 
can be shown to be determinate in the same manner. For 

A = 0 the problem becomes indeterminate, that is the induced 
field cannot in general be specified in terms of the 


inducing field, as discussed by Price (1950). 


2.4 The Uniform Source Two-Dimensional Solution 
The form of the total electric or magnetic field 
for the uniform source for either the E or H polarization 


case can be found from the solution of equation (2.57) or 


(2.58). Considering first the H-polarization case, fora 
oH 

uniform conductor ae = 0 and (2.58) reduces to 
aebilee Jes 
ee = iwy oH, for z>0 (2..85) 


This, together with the boundary condition that He tends to 
zero aS z + @ leads to the solution 
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inside the conductor. 

By virtue of equations (2.52) and (2.53), i in 
the region outside the conductor, whereo=0, iS independent 
ef both y and z and is, therefore, uniform throughout this 
region. The H field is continuous across the surface boun- 
dary, and so ie is constant everywhere just inside the 


surface of the conductor. This implies that oH / oy is zero 
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just inside this surface and therefore from (2.53) that ob 
is zero. This agrees with the fact that normal component of 
Curnent must be zero at the surface.) /#1f (2.58) 4s solved 
with the appropriate value of © and the boundary conditions 
on Hs then the required conditions on E in the conductor 
will be satisfied. 

In the E polarization case, only the field 


components EL: Ho and ue are involved, and EY Satisfies 


y JE 
(225/). Again, for a uniform conductor, ATE =“ Oramdeie a7) 
reduces to 
9°E 
> = LOMO a iar (2587) 
dZ 
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E ure Thirve 


fot 7-0 (268') 


From equation (2.56) it can be seen that H, is zero. By 
substituting (2.88) into (2.55) the solution for HY can be 


found inside the conductor. 
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In the non-conducting region oH a2 goes to zero by virtue 


of equation (2.51). This means that HY must be constant 
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(say Ho) with respect to z outside the conductor. Since HY, 
is continuous across the surface, from (2.89) we have 


(wn,0) 72 


1 
pe” Og ame ae ee (2.90) 


From (2.55) in the non-conductor 


aE 
eg ar TwOH OH, V2.0") 
eid eee ee aa ey (2.92) 


Within the conductor all field quantities vanish as z+, 


Other boundary conditions are that EY. H and He are 


ah 
Pomenuous ate Z=0.04ASs Ze -.°,, Ele goes, to infinity, 1dbis 
result is in accordance with an energy source which must 
exist at or beyond the plane z = -|h| . This analysis, 
chi confined to uniform subsurfaces, can be extended to 
multi-layered cases. 

In what follows, the E-polarization case is used 
throughout for the uniform source boundary conditions. By 
solving for the electric field in the two-dimensional 
ealculations. these electric. field solutions may be used 


directly in the Lines and Jones (1973) finite differences 


method to provide boundary and initial conditions for the 


mesh. 
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2.9 The Two-Dimensional Non-Uniform Source Solution 

In this section a two-dimensional, non-uniform 
Sheet current source with Gaussian current intensity distri- 
bution is discussed, and the solution for the electromag- 
netic field over a layered conductor is derived. In the 
next section (2.6) a method of using this solution to 
generate any arbitrary two-dimensional source configuration 
Psi ven. 

The coordinate system used is shown in Fig. 2.2. 
The layered conductor occupies the half-space z > 0 as 
Shown. The sheet current source is assumed to flow parallel 
to the surface of the earth at a height z = -|h|. Following 
Peltier and Hermance (1971) and from equation (2.53), it can 
be shown that for a horizontally layered earth and for a 
symmetric current source the elementary electric field in 


the n-th layer is of the form 


~Oy2 NES 
Ee = {A, e Pape. sh ete cee (2.93) 
Lely : 
where ox =" + iw, Fp 


Bid Geng pe. and o, are respectively the spatial wave number, 
f°) 
the angular frequency, the permeability of free space and 


the conductivity of the n-th layer. For a sheet current with 


a Gaussian distribution function 
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The coordinate system, source position and 
layered model. 
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54 
where Ne is the maximum intensity and « is the standard 
deviation of the source. 

The quantities Ae and Bi are determined by the 
conductivity configuration and the solution for the electric 
field is completed by taking the sum of Ey over all A: 


=028 Z GZ 


Pye) { fA en”) +B e " } Cfabeostiyjan . 202.959 
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After calculating EL the magnetic components may 


be determined: 
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(2.96) 
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Pely sz) = (ivyou,) if {AFA e + Be J CIA)! sinfaydad) . 
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The half-space above the Earth is considered as 
mene first layer, and if A, = ] then C(A) is a frequency 
dependent parameter of the source. B, may be thought of as 
a reflection coefficient which represents the contribution 
to the total electric field at the surface by the subsurface 
layers. For a three layered conductor By =.0 so «that EY 
will not become infinite with depth. The other A, and BY 


values may be determined Since the tangential components of 


Ee" and H" must be continuous across layer boundaries. This 
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where qd, and d., are the thicknesses of the first and second 
conducting layers respectively. Also, we have 
(1 + B,)Q 
A, = ——-— ; 


(2.99) 


and By = Ay Q, 
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29 
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A, = tA. + Bae 
This general method may be extended to more layers 
if desired. It should also be pointed out that the conver- 
gence of the field integrals is discussed by Hermance and 
Peltier (1970) and Peltier and Hermance (1971). 
The above solutions are limited in that only source 


configurations for which the source coefficients C(A) as 
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described by Peltier and Hermance (1971) are obtained 
analytically could be modelled. It is clear that more 
General source configurations exist, particularly in: high 
latitude regions (Kisabeth and Rostoker, 1971 ; 


Kisabeth, 1972). 


2.6 A Method for Obtaining an Arbitrary Two-Dimensional 


source 
A method in which any arbitrary two-dimensional 

Source configuration may be represented is desirable and in 
this section such a method is given. This method provides 
field values for a layered earth, and more general problems 
may be studied by using these values as boundary conditions 
for a numerical technique. The computer program listing for 
this method and description of its use are given in Appendix 
A. 

A general source intensity distribution may be 
constructed from a number of elementary sources by super- 
position. A Gaussian intensity distribution of small half- 
width (1 - 10 km) has been chosen here as the elementary 
Source, since C(A) can be obtained analytically for this 
distribution. A number of these elementary Gaussian sources 
of equal maximum intensity are then spatially shifted so 
that a rectangular current intensity distribution may be 
approximated. The electric and magnetic field values are 


then obtained by superimposing the spatially shifted fields 


of the elemental sources. 
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In this way an arbitrary current intensity distri- 
bution may be approximated in a piecewise continuous manner 
by weighting a number of spatially shifted rectangular 
current distributions in accordance with the trapezoidal 
rule where the field components are spatially shifted and 
Superimposed. The solution for a two-dimensional source of 
arbitrary intensity distribution over a layered Earth may 
thus be approximated. This is equivalent to convolving the 
arbitrary current distribution with the Green's Function of 
the rectangular current distribution. 

| The arbitrary current source could be approximated 
by directly shifting and summing the elemental Gaussian 
field solutions using a criterionsuch as a least squares 
process to weight the Gaussians. However, this would increase 
the number of computations required by a factor equal to the 
number of Gaussians used to construct the rectangular current 
source distribution and so 1s not desirable. 

Because of its simplicity, the method is readily 
applied to modelling problems. The width of the rectangular 
distribution and therefore the half-width of the elemental 
Gaussian can be easily estimated by linearly approximating 
the arbitrary current source to the desired accuracy. AS 
many elemental Gaussians as desired can be used to represent 


the rectangular distribution and for the examples here 3 and 


31 Gaussians have been used. 
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2.6.1 Comparison of Piecewise Continuous Representa- 
tions and Analytically Integrated Kernal 
Calculations 


In order to obtain an indication of how well the 
elemental Gaussian and rectangular representation of a 
source will approximate the solution for a case in which 
C(A) may be represented analytically several comparisons are 
made. It is advantageous to represent the rectangular cur- 
rent distribution using as few elemental Gaussian currents 
as possible. To determine the number that may be used a com- 
parison of the field values calculated by using 3 elemental 
Gaussians with those calculated by 31 elemental Gaussians to 
represent the same rectangular source was made. Calculations 
of the field values at the surface of a two-layered Earth in 
which a 50 km layer of 100 ohm-m resistivity overlays a 
half-space of resistivity 10 ohm-m at three frequencies 
(1.0 Hz, 0.01 Hz and 0.0001 Hz) were made for Gaussian sheet 
current intensity distributions with half-widths equal to 
0.968 km and 10.0 km. A rectangular current distribution of 
30 km width was approximated by 31 of the 0.968 km elemen- 
tal Gaussian currents spatially shifted one kilometer and 
compared with a representation of the 30 km rectangular 
current composed of 3 elemental Gaussian currents of 10 km 
half-width shifted one half-width. The results are shown 
in Tables 2.1 and 2.2. The horizontal components are 
normalized at the origin in both tables, whereas H, is 


normalized at 15 km in Table 2.1 and 400 km in Table 2.2. 
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Comparison of Gaussian representations of rectangular sources near the origin 


Distance 

from 
origin 
(km) 

0 

15 

30 

60 


Distance 
from 
origin 
(km) 
0 
15 
~ 30 
60 


Distance 
from 
origin 
(km) 
0 
15 
30 


Ex 


0-9833 
09-360 
0-7839 
0-6158 
0-4729 
0:3637 


31 
1:0 
0:9891 
0-9575 
0-8498 
0-7168 
0-5882 
0:4778 


Frequency: 1 Hz 


Hy 


31 
1-0 
0-9833 
0:9359 
0: 7837 
‘)-6156 
0:4726 
0:3635 


Frequency: 0-01 Hz 


31 
1:0 
0:9868 
0:9492 
0: 8258 
0: 6833 
0:5539 
0:4477 


Frequency: 0-0001 Hz 


31 
1:0 
0-9939 
0-9760 
0:9131 
0-8301 
0:7427 
0:°6595 


A, 
31 


-1279 x 1075 


A, 
31 


4178 x 10-° 
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Comparison of Gaussian representations of’ rectangular sources at great distances 
from the origin 


Frequency: t Hz 


Distance 
from 
origin E. H, A, 
(km) 3 31 3 3t 3 31 
0 1-0 1:0 1-0 1:0 0:0 0° 8233 x 1075 
400 0-07232 0: 7293 0:-07204 0:07265 1°0 1:0 
800 0:01829 0-01842 0-01796 0-01809 0- 1391 0-1391 
1200 0:-007588 0-007640 0:007248  0-007297 0:04297 0-04274 
Frequency: 0-01 Hz 
Distance 
from 
Origin E, H, A, 
(km) 3 31 3 3 2 31 
0 1-0 1:0 1-0 1:0 0:0 0:1222x 1075 
400 0:1110 0-1116 0- 1052 0-1059 1-0 1-0 
800 0-02912 0-02926 0:02742 0:02759 0: 1472 0-1470 
1200 0:01240 0:01246 0:01143 0:01149 0:04514 0°04504 
Frequency: 0-0001 Hz 
Distance 
fror 
origi E, HA, A, 
(km) 3 31 3 31 3 31 
9 1-0 1-0 1-0 1-0 0-0 0-1133 x 10-6 
409 0: 2493 0:2501 0-1707 0-1718 1-0 1:0 
800 0:-07675 0-07697 0:05785 0-05820 0:2019 0-2018 
1200 0:03455 0:03465 0:02609 0:02625 0:06467 0-06463 
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From the tables it is clear that the three 10 km elemental 
Gaussian currents well represent the rectangular current 
Bretri bution In the calculation, of the field values. 

Also, in Tables 2.3 and 2.4 the fields associated 
with a Gaussian sheet current of 240 km half-width calculated 
in the manner of Peltier and Hermance (1971) and Hibbs and 
Jones (1973a) are compared with those of a 49 element 
piecewise continuous source constructed from 30 km rectangu- 
lar current sources employing the two approximations as used 
above. The same subsurface and frequencies as in Tables 2.1 
and 2.2 were used. EY and Hy are normalized at the origin 
in both tables while H, is normalized at 15 km in Table 2.3 
and 400 km in Table 2.4. 

Table 2.3 shows that for distances of 150 km or 
less from the origin the two solutions are nearly identical. 
Table 2.4 indicates that the two solutions differ as the 
distance from the origin increases. This is because the 
piecewise continuous source is terminated at /20 km whereas 
the Gausswanesource continues: *toginfinity. gahisseas particu- 
larly evident in the H, component which one would expect to 
increase when the current intensity changes rapidly. 

A further comparison with previous results was 
made for a sheet current source of intensity distribution 
(y-b) igi! u(y-b) where b = - 480 km and u(y-b) is 
the unit step function. Using the present method, the 


source was approximated by forty-nine rectangular sources 
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of 30 km width. Each rectangular source was approximated by 
three elemental Gaussians of half-width 10 km and Spatially 
shifted one half-width. These calculations were made over a 
range of periods from 1.0 sec to 10° sec for the same sub- 
Surface as above. The apparent resistivity (calculated in 
the manner of Cagniard, 1953) and phase curves for this case 
are shown in Figs. 2.3 and 2.4. The apparent resistivity 
curves are calculated for points at intervals of 200 km 


3 to 10° 


over the range from - 10 km. The phase calculations 
illustrate the difference between the phase of EY and the 
phase of HY, for the various frequencies at the different 
positions on the surface of the Earth beneath the source. 


3 and 10° sec 


The field profiles of E> Hy, and ie FOF 1 Oe 0 
period are shown in Fig. 2.5. 

The apparent resistivity curves of Figs. 2.3 and 
2.4 illustrate the typical spreading as a function of 
position at longer periods and compare well with those 
obtained previously (Hibbs and Jones, 1973b). A.comparison 
of apparent resistivity curves is a good check on the 
accuracy of the method. Highly accurate field values are 
required to calculate the apparent resistivity curves since 
the field values are squared and a ratio taken. If the 
field values are inaccurate this produces spreading of the 
apparent resistivity curves at short period as well as 


incorrect values at long period. Also the field profiles 


for E. and H. of Fig. 2.5 compare favourably with previous 
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Figure 2.3 


Apparent resistivity and phase curves for th 
(y-b) fe[p(y-b) Ju(y-b) source over a two- 
layered conductor with a 50 km thick upper 1] 
of resistivity 100 ohm-m overlying a half-sp 
of resistivity 10 ohm-m. The curves are for 
the negative y region and are calculated for 
points along the surface beneath the source. 
Square >: -1000 km 
circle > - 800 km 
triangle : - 600 km 
3 > - 400 km 
X > - 200 km 


diamond : at origin 
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Figure 2.4 As Figure 2.3 but for positive y region. 
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Field profiles along the surface of the earth 
for three periods for model as in Figure 2.3. | 


The field components are normalized with respec 


to the point at the origin. 
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work. However, the vertical magnetic component (Ho) 
increases near 10° km. This is because the current source 
is terminated at + 1010 km whereas previously the source 
continued to infinity. The results shown in these plots, 

as well as the apparent resistivity ones of Figs. 2.3.and 
2.4 reflect the non-symmetric character of the source. 

Fig. 2.6 gives the apparent resistivity and phase 
curves for the elemental Gaussian itself. The Gaussian 
source is symmetric and the results for + y and - y are 
identical. The spreading of the curves at longer period as 
a function of position is again evident. At short periods 
the phase difference between EY and Hy is 45 degrees at all 
positions beneath the source. However, as the period 
increases this phase difference becomes dependent on the 
position relative to the source. This effect is evident 
for all the non-uniform source models. The field component 
profiles are shown in Fig. 2.7. 

The second conductivity configuration considered 
was that of a two-layered Earth with upper layer 4 km thick 
and resistivity 0.25 ohm-m with an underlying half-space of 
resistivity 250 ohm-m. The source was constructed in the 
Same manner as in the previous model. The apparent resistivi- 
ty and phase curves are shown in Figs. 2.8 and 2.9 and the 
field component profiles in Fig. 2.10. The apparent 


resistivity curves again compare well with those obtained 


previously by Hibbs and Jones (1973b). 
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Figure 2.6 Apparent resistivity and phase curves for the 
elemental Gaussian source over a two-layered 
conductor with 50 km thick upper layer of 
resistivity 100 ohm-m overlaying a half-space 
resistivity 10 ohm-m. The curves are for the 
positive y region and are calculated for point: 


along the surface beneath the source. 
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Figure 2.7 Field profiles for model as in Figure 2.6. a | 


and Hy are normalized with respect to the origi! 


and H. with respect to the point at y = 200 km. 
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Figure 2.8 


Apparent resistivity and phase curves for the 
(y-b) e r(y-b) u(y-b) source over a two Tayere 
conductor with 5 km thick upper layer of resist 
tivity 0.25 ohm-m and half-space of 250 ohm-m.- 
These curves are for the negative y region and 
are calculated for points along the surface 


beneath the source. 
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Figure 2.9 As Figure 2.8, but for positive y region. 
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Figure 2.10 


Field profiles for model as in Figure 2.8. 
All components are normalized with respect to. 
the origin. The results for the 10 sec and th 
10° sec curves are nearly the same, and these 
two curves overlap. 
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2.6.2 An Arbitrary Non-Integrable Source Configuration 

A third more general model was considered. The 
source for this model was composed of fourteen rectangles, 
each 30 km in width and with intensity coefficients as given 
in Table 2.5. Each of the rectangles was approximated by 
three elemental Gaussians of 10 km half-width. The conduc- 
tivity configuration was that of a two-layered Earth with 
upper layer 50 km deep and resistivity 100 ohm-m. The 
second layer was a half-space of resistivity 10 ohm-m. The 
apparent resistivity and phase curves are given in Figs. 
eal. wand 2.12 and the figldncomponent profiles in Fig. 2.13. 

For this third arbitrary model the spreading of the 
apparent resistivity curves for different positions underthe 
source at long periods is evident as in the other two models. 

By superimposing elemental Gaussian solutions to 
construct rectangular current segments any general source 
configuration can be approximated. This provides much 
flexibility in the construction of sources to aid in the 
removal of source effects from array data and also in the 
study of source configurations from ground measurements when 
the conductivity of the earth is considered. This approach 
may also be used to provide boundary conditions for more 
general laterally inhomogeneous earth model calculations. 
A description of the actual program and a discussion of its 


use is given in Appendix A. 
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Table 2.5 Source intensity coefficients for the arbitrary 


source. — 
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Figure 2.11 


Apparent resistivity and phase curves for the 


arbitrary source over a two-layered conductor . 


with a 50 km thick upper layer of resistivity 


100 ohm-m overlaying a half-space of resistivi. 


10 ohm-m. The curves are for the negative y 


region and are calculated for points along the 


surface beneath the source. 
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Figure 2.12 As Fig. 2.11; but for the positive y region. 
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Figure 2.13 Field profiles along the surface of the earth 
| 


for three period? for model of Fig. 2am 
field components are normalized with respect to 
the point at the origin. 
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2./ The Horizontal Magnetic Dipole Source Solution for a 

The horizontal magnetic dipole solution presented 
here was originally formulated by Grant and West (1965) for a 
magnetic dipole above a homogeneous earth. The solution for 
the polar magnetic fields was derived for a monopole and 
then by simple differentiation the solution for a magnetic 
dipole was obtained. Ramaswamy (1973) extended this work to 
include the effect of a two layered conducting earth. 

Consider a semi-infinite, tnree-dimensional, two 
layered conducting earth occupying the half space z >0O of 
the rectangular coordinate system shown in Fig. 2.14. The 
first layer is of conductivity o, and depth d, and the second 
layer is semi-infinite with conductivity Oo: The upper half- 
Space is a free space region. The permeability of both the 
conductors and the upper half-space is equal to that of free 
Space. From the mathematical model shown in Fig. 2.14, a 
magnetic pole situated at (0,0,-h) has cylindrical symmetry 
and cylindrical coordinates (r,®,z) are used. These are 


related to the Cartesian coordinates (x,y,z) by 
4 


x = r cos 0 and: y = r sin 6 : (2.100) 


The vector potential, A, for the magnetic pole will have 


only a @ component due to symmetry which will be a function 
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Figure 2.14 The magnetic dipole and the coordinate system. 
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In the three conductivity regions considered, the 
source currents are excluded so that equation (2.10) must be 
sarustied in each ‘région. “(r,z) “automatically satisfies 


mreunelation Vv + =-0Y “ Thevetore, 


poten oe - in wo,m(r sz) = 0 27 0a 


where Wo? and o, are the permeability of free space, the 


n 
angular frequency and the conductivity of the n-th layer 
respectively. The time dependence is sinusoidal and the 
factor anal is understood in all field quantities. 
Solutions of (2.101) are of the form 
[+ (6° + injwo,)!/?2] 

e Ji (rf). G22 i029 

where J,(r.é) is a Bessel function of the first Kind ‘of 


order 1, and 
Keene ts jie Os eae 

2 Yi (r B} (2e810'3)) 
where V1 (16) is the Bessel function of the second kind of 
order 1. Since the fields must be finite for r=0, z#0,only 
solutions of the form (2.102) are permitted. The plus and 
minus sign in the exponential are chosen so as to reject 
solutions which are not finite as z > + ~. The complete sol- 
ution of (2.101) is obtained by superimposing all possible 
solutions of the type (2.102). The magnetic field of the 


pole is obtained from relation (20er7 ) 
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DOT ever atc?) a (21.00) 


where 8 is the uhtt vector in the. 6 direction. 

The horizontal magnetic dipole solutionis obtained 
by considering magnetic poles of strength (+p) and (-p) vary- 
ing sinusoidally with time whicn are situated at (X,,Y) .-h) 
and (XxX), Y, + AX, -h) respectively. The total field at any 
observational point (x,y,z) is obtained by superimposing the 
fields of the individual poles. In the limit as AX +> 0 the 
total field will reduce to that of a dipole at (X, > ie -h) 
with its dipole moment directed in the negative y-direction, 
provided that AX approaches zero in such a way that p 
remains constant. That,is, if URE and ae Sener 
the fields at a point (x,y,z) due to poles +p and -p, 
respectively, then the field at (x,y,z) can be written as 
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The dipole moment m is given by tne relation 
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The solution for I depends only on the horizontal distance 


between the pole and the field point, therefore ar = - a : 
] 
And 
M or 
he ee ae 2.109 
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Wye treid, J, is related to thé vector potential by (2.17) 
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Using the elementary coordinate transformation 


um = r cos 9 and ay = i Sin 8 the x and y components of 


the dipole field are 
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The solution for the electric field is found by taking the 


curd of Hh. 
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The solutions for the magnetic and electric fields 


in the semi-infinite layer are: 
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The magnetic field solutions for the free space 


region are found in the same manner as before. 
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In the free space region the electric field cannot be ob- 
tained through the use of the curl equation since o = 0 


there. However, the relation 
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These equations must be integrated in order to find the 
ervectric field. Since ES vanishes inside the conductor, the 
conductor acts like a perfect reflector (infinite conduc- 
tivity) so far as E, is concerned. Thus geometric imaging 


of the dipole can be used to find the solution for ES: 
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substituting (2.136) into’ (2.133) and integrating the 


resulting equation with respect to z gives 
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The constant of integration is zero since Ey in 62-132).must 
match Ey given by (2.118) in the region 0 < z < d at the 
boundary z=0. Similarly substituting into (2.134) and 


integrating with respect to z 
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In this way, the fields associated with a horizon- 
tal dipole above a layered conductor may be determined. The 
foregoing solutions were used for all boundary conditions 


involving horizontal magnetic dipoles in this work. 


2.8 Three-Dimensional Finite Difference Numerical Technique 
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dimensional methods. However, some conductivity anomalies 
Such as the Alert, Tucson and Japanese anomalies require a 
three-dimensional treatment. Also any embedded anomaly with 
approximately equal spatial dimensions requires a three- 
dimensional treatment. 

Jones and Pascoe(1972) treated the local induction 
problem in three dimensions for an imbedded cube. Maxwell's 
equations were solved over a three dimensional uniform grid 
by the use of finite difference approximations and the Gauss- 
Seidel iterative method. Lines and Jones (1973) extended 
the work to a grid with variable dimensions so that more 
general anomalies could be accommodated and so that the 
external boundary conditions could be better satisfied. 

The basic method is described in Lines (1972) and is 
Summarized here. 

If the curl of equation (2.2) is taken and 

substituted into equation (2.1) then the resulting equation 
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mearee space. All Components of Eowary With oxet iv. cand. -z 
in Cartesian coordinates. 
The above vector equation in E can be rewritten 


as three scalar equations 
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These three equations are solved Simul taneously FO? E> EY 
and ES 

In the discussion of the numerical method, only 
(2.140) will be considered. A similar procedure is used in 
evaluating (2.141) and (2.142). The three-dimensional finite 
difference grid is shown in Fig. 2.15. The numerical formu- 
jacion forepaint '0' is considered. The distance to adjacent 
points in the grid are in general unequal. Eight conduc- 
tivity regions: surround the point, “0°. jhe subscripts of 
the electric field components refer to the points on the 
grid ‘of “Fig. 2.75 ‘at which ethe components, ane evaluated. 

In (2.140), the finite difference expressions for 
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Figure 2.15 Finite difference grid (from Lines, 1972). 
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2. 144) 


where the central differences have been used for (dE /ay). 
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Similarly for (2.141) and (2.142) 
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dimensional work of Jones and Pascoe (1971), the 'fictitious' 
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the internal boundary conditions. The condition that the 
tangential components of E and H must be continuous eliminate 
2/3 of the ‘fictitious values'. A problem arises when E is 
normal to a boundary for which o jis discontinuous. J must 
be continuous across the boundary. If o is discontinuous 
at a boundary, then E must also be discontinuous. However, 
this cannot be represented by a point value at the boundary. 
Therefore, the average of E on each side of the boundary is 
used. When expressions of (2.149), (2.150) and (2.151) are 
Summed over all eight regions surrounding the point 'O' and 
the boundary conditions are applied to E, the resulting 
expressions are of the same form as in (2.149), (2.150) 
and (2.151) with o being replaced by o the average of all 
the Pr iniicii os of the surrounding regions. This implies 
that discontinuities in conductivity are represented by tran- 
sition zones between the regions of different conductivity. 
The resulting finite difference equations (2.149), 
(2.150) and (2.151) are solved by using the Gauss-Seidel 
relaxation technique for simultaneous linear equations. 
Initial values are required throughout the interior 
ef the grid,.to, stant, the jterative: process. ..As. shown jim the 
preceding sections of this chapter, these initial values may 
be determined analytically for a layered half-space. In the 
next section various source configurations are considered 


and the boundary conditions for the finite difference mesh 
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are determined. The solution of the local three-dimensional 
problem is then obtained for various embedded cConductLvity 


anomalies and frequencies. 


9.9 Electromagnetic Induction in Three-Dimensional 


Structures for Various Source Fields 

In this section the numerical method is used to 
calculate the electromagnetic fields associated with three- 
dimensional conductivity anomalies. Three source field- 
conductivity structure combinations are considered. A 
uniform source above an island structure as well as a two- 
dimensional non-uniform source above an L-shaped embedded 
Structure and a horizontal dipole above an embedded square 
anomaly are studied. Profiles of selected electromagnetic 
field quantities and ratios are presented along with the 
appropriate phases for various positions relative to the 
anomalous structure. Three dimensional figures of the 
amplitudes of the electric and magnetic field components 
are given for each source field-conductivity structure 
combination. The source field effect is clearly evident in 
‘mall profiles. 

The electromagnetic fields due to one, two, and 
three-dimensional sources over layered laterally uniform 
conductors are calculated as discussed previously. For the 
One-dimensional source the analytical solution for the 


layered earth as in section 2.4, is used here. For the 
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two-dimensional source configuration the Fourier series 
method as described in section 2.6 is used, and for the 
three-dimensional dipole source the method of Ramaswamy 
(1973) as described in section 2.7 is used. As pointed out 
before, the electric field values thus calculated are used 
to provide boundary and initial conditions for the finite 
difference mesh of Section 2.8 and from this the 

electric fields due to three-dimensional embedded conduc- 
tivity anomalies are determined. The magnetic field 
quantities are then obtained from the electric field solu- 


tions by using the appropriate finite difference equation. 


ao One-Dimensional Source Field 

The one-dimensional source field consists of a 
uniform sheet current flowing above a layered conductor. 
In the present work, the E-polarization case is considered 
and the electric field component is polarized in the x- 
direction with the magnetic field component polarized in 
the y-direction. 

In the model considered, a uniform sheet current 
flows above a conductivity inhomogeneity embedded in a two- 
layered semi-infinite conducting region with a plane boun- 
dary, which approximates an island in the deep ocean for 
geomagnetic bay type disturbances. The first layer is 4 km 
in depth and of conductivity 4 mho/m. The second layer 
extends to infinity and has conductivity 0.001 mho/m. 


The embedded anomaly is a 6 km square’ intrusion 
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of the lower conductive layer which penetrates the top layer 
to the surface of the conductor. The calculation was done 
for a frequency of 0.001 Hz. 

From the foregoing solution the profiles of Ee 
Hy Ho HO fH, and the apparent resistivity (calculated in 
the manner of Cagniard, 1953), as well as their respective 
phases were calculated for 3 profiles at different positions 
relative to the anomaly as in Fig. 2.16. The profiles were 
calculated over the range - 58 km to + 58 km in the y 
direction. Each field component or ratio with the exception 
of the apparent resistivity and }H O/H | was normalized with 


respect to its value at -58 km. All phases are normalized 


to zero at - 58 km. 


The first set of profiles is taken 8 km from the 
center of the anomaly in the negative x direction and is 
given in Fig. 2.17. The normalized electric field amplitude 
Es exhibits a slight depression in the field due to the 
concentration of charge on the boundary of the anomaly. The 
normalized horizontal magnetic field [HA | is reasonably 
constant across this profile and changes by only about 4 
percent. Normalized |H,| and therefore JH O/H, | 
vary in the same manner over the anomaly due to the constant 
nature of HY Since both rey and [H | are nearly constant 
the apparent resistivity (o,) curve is nearly constant at a 
value of 1 ohm-m. The value differs from the actual resis- 


tivity of .25 ohm-m because of the long period used since 
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Figure 2.16 Anomalous structures with profile positions 
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profiles with normalized phases for one- 
dimensional source field and profile 1. 
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the calculated apparent resistivity value is influenced by 
the lower conductivity region. The phases are nearly 
constant, varying by only a few degrees, for all profiles 
except ie, and ol ty As would be expected, the phase of HE 
and therefore HOfH shifts by 180° over the center of the 
anomaly, which is consistent with the change of sign of ue 
there. 

The second set of profiles, taken along the line 
XY = see2 KieetSugiven in Fig. 2.18 and cleanly .illustrates the 
effects of the anomaly. JE, | increases over the anomaly and 
JH decreases. |H, | and JH O/H | markedly increase over the 
anomaly. The apparent resistivity now shows a small 
increase in value over the anomaly and its value is two 
orders of magnitude below the actual resistivity value. The 
oe of the anomaly is insufficient to cause an appreciable 
change in the apparent resistivity. The phase of the elec- 
tric field, EY: changes by a maximum of 23 degrees over the 
anomaly when compared to the phase of the electric field at 
-58 km. The phase of the apparent resistivity is similar 
to the electric field phase since the phase of Hy is nearly 
constant. Again the phases of He and H O/H change by 180 
degrees as the center of the anomaly is crossed. 

The third set of profiles, taken across the center 
of the tanoma ly jis given,in Fig. 2.194. This set 4s almost 
identical with the second set of profiles except that eo 


is more sharply peaked as would be expected. 
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Three-dimensional amplitude profiles of the six 
electromagnetic field components are presented in P¥Gn “2.20 < 
The decrease in the re component off the anomaly in the 
x-direction due to the concentration of charge at the boun- 
dary of the anomaly is clearly evident. Over the anomaly 
|E,| increases while |H | decreases. [JH | increases over 
the edges“ot they ahomaly in the + Y and ’-/ Widirectionsx 
Near the corners of the anomaly both [EY | and ga increase 
as would be expected due to bending of the fields in these 
reqvons -//)The qvertical electric field, [E,|» increases on 
both sides of the anomaly in the + X and - X directions. 
Wis is because the current is deflected vertically by the 


island structure. 


Dehie Two-Dimensional Source Field 

The two-dimensional source field is obtained from 
a non-uniform sheet current flowing above a layered conduc- 
tor. Again, in the present work, only the E-polarization 
case is considered since it is most easily adapted to the 
formulation of the solution for the electromagnetic fields 
of a non-uniform sheet current source. For this case both 
the horizontal magnetic field component, Hs and the verti- 
cal magnetic field component, H,, will be present in the 
solution for a layered conducting subsurface. 

In this work, a sheet current of rectangular 
intensity distribution, simulating an ionospheric 


electrojet, which is 1440 km in width and centered 
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840 km from the center of the L-shaped embedded anomaly 
(Pra: F2°V6Eby "Ts \Gonstdered. /This sheet current flows at a 
height of 110 km above the surface of the conductor. The 
conductivity anomaly is embedded in a two-layered, semi- 
infinite conducting region with a plane boundary which is 
chosen to represent an average continental vertical conduc- 
tivity distribution. The first layer is 50 km in thickness 
and of conductivity 0.01 mho/m. The second layer extends 
to infinity and has a conductivity of 0.1 mho/m. The 
embedded anomaly is L-Shaped, as in Fig. 2.16b, and has a 
depth of 3 km. The conductivity of the anomaly is 0.0001 


mho/m and the calculation is done for a frequency of 0.1 Hz 
in order to enhance the magnitude of the anomalous field. 


The boundary conditions for the finite difference 
mesh were determined by using the general method for two- 
dimensional sheet current distributions over layered conduc- 
tors as described in section 2.6 above. These values were 
then substituted into the numerical technique of Section 
2.8 and the electromagnetic field solutions obtained. 

The three profile positions chosen for this anoma- 
ly are shown in Fig. 2.16b. The profiles are calculated 
over the range - 38 km to + 38 km in the y-direction. Each 
field component or ratio with the exception of the apparent 
resistivity and [H/H, | is normalized with respect to its 
value at -38 km. All phases are normalized to zero at - 38 km. 

The first set of profiles, which is taken 6 km 


from the center of the anomaly in the X-direction is shown 
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in Fig. 2.21. The sloping character of the source field is 
clearly evident in the (eet EY | and [RiP \pnetites.. Since 
ime-conductivity~contrast is oly 100 te 1, the effect of 
the anomaly at this profile position is not evident for Hey 
and |H| at this frequency. Also, as will be seen later, 
the rapid change in the source field tends to screen the 
effect of the anomaly when the fields are normalized. The 
effect of the anomaly can be clearly seen in the |H | and 
JH /H | profriles. The apparent resistivity is constant for 
this profile and corresponds in value to the conductivity of 
the upper layer. The phases are sensibly uniform for all 
quantities for this profile position although the expanded 
scale of the diagram amplifies any change in phase. 

The second set of ‘profiles (Fig. 2.22), A= taken 
2 km from the center of the anomaly in the X-direction. 
These show effects due to the anomaly. The anomaly is just 
evident in the Ry profile while in the |H | profile the 
dipping of the profile due to the anomaly is apparent. The 
JH O/H | profiles clearly show the effect of the anomalous 
structure. As before, the apparent resistivity profile 
snows a slight increase due to the anomaly. However, the 
size of the structure is too small for the apparent resis- 
tivity to have a value equal to the resistivity of the 
anomaly which is 10000. ohm-m. A noticeable change in the 


phases of both H, and H/Hy is now evident over the anomaly 
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while the other phases remain relatively constant. 

The third set of profiles, Given in/Fig. 2.28 is 
much the same in form as the second set of profiles. However, 
the effect of the anomaly is apparent over a greater distance 
than before due to the increase in the width of the anomaly 
atithis profile position. Alsa@;} more clearly evident in the 
JH | and JH 7H | profiles is the induction effect of the H, 
component of the source field. The peaks in the [Ho profile 
are not symmetric and differ in character. 

In Fig. 2.24 the three-dimensional amplitude 
profiles of the six electromagnetic field components Ae 
presented. The figures are distorted near the edges since 
a non-uniform grid was used in calculating the field solu- 
tions while a uniform grid is used in the three-dimensional 
presentation. The anomalous electric field in the x-direc- 
tT On ; JEL > is barely discernible due to the rapid spatial 
change of the source field. All the other electromagnetic 
field components clearly show the effects of the non- 
symmetric-structure.- The corners -of-the-structures—are 
identified with large increases in [EY and [HI . HI 
dips over_the-structure and the;edges of the structure 


are outlined by et and eae 


2.9.3 Three-Dimensional Source Field 
A three-dimensional source field can be obtained 
by using a horizontal magnetic dipole as a source as out- 


lined in section 2./7 above. 
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In the model considered a horizontal magnetic di- 
pole was placed 17 km in the negative y-direction from the 
center of the anomaly (see Fig.2.16c) at a height of 150 km 
above the surface of the anomaly. The dipole is polarized 
in the negative y-direction. The conductivity anomaly is 
embedded in a two-layered, semi-infinite conducting region 
with a plane boundary. The first layer is 128 km in depth 


and of conductivity 0.21 x 107° 


mho/m, The second layer is 
of infinite depth and has conductivity 0.8 mho/m. The 
embedded anomaly (see Fig. 2.16c) is 16 km square and 4 km 
deep with its top at the surface of the conducting region 
and has the same conductivity as the lower layer. This 
calculation was carried out for a frequency of 0.075 Hz. 

The two profile positions for this anomaly are 
shown in Fig. 2.16c. The profiles are calculated over the 
range y=-42 km to y=+ 42 km. Each field component or ratio 
with the exception of the apparent resistivity and |H,/H)| is 
normalized with respect to its value at - 42 km. All phases — 
are normalized to zero at - 42 km. 

The first set of profiles which is taken 24 km 
from the center of the anomaly in the X-direction is shown 
in Fig. 2.25. The non-uniformity of the source field is 
clearly evident in all the amplitude profiles. For this 
frequency and conductivity contrast the anomalous structure 
is evident at this distance away from the anomaly. 
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profiles with normalized phases for three- 
dimensional source field and profile 1. 


— Total field due to anomaly 
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in the fe profile and a depression in the [HI 
profile. The boundaries of the anomaly are shown in the 

[HS | and JH O/H | profile by slight peaks. The apparent 
resistivity profile exhibits a slight increase, opposite to 
what would be expected, because of the increased electric 
field in that region caused by the charge concentration. 
The phases are nearly constant with the exception of the Ho 
and Hf phases which go through a change of 180 degrees 
as the vertical component of the source field goes through 
zero. 

The second set of profiles, taken across the center 
of the anomaly, is shown in Fig. 2.26. In these profiles 
the effect of the anomaly is clearly displayed in all 
amplitude components and their phases. A marked depression 
in the RES | profile is evident. JH | increases rapidly over 
the anomaly and the }H | and JH 7H | profiles exhibit 
two peaks over the edges of the anomaly. The apparent resis- 
tivity varies from approximately the value of the resistivi- 
ty of the first layer to the value of the resistivity of the 


anomaly. Increases in phase are noted for Es H and Pa 


y. 
while the phases of H, and nih change by 180° as the source 


field goes through zero. 
In Fig. 2.27 the three-dimensional amplitude 
profiles of the six electromagnetic field components are 


presented for the three-dimensional source field case. A 
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large depression in Boe is present over the anomaly while 
just outside the anomaly the field increases slightly due to 
the varying surface charge on the boundary between the 
conductivity regions. }H.| increases over the anomaly and 
the corners of the anomaly are marked by increases in |H, | 
and [EVI The sides of the anomaly are defined by FES and 
[Ho I. The peaks of de. are not of equal height because of 


the non-uniform source field. 
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A COMPARISON OF ELECTROMAGNETIC ANALOGUE MODEL 
MEASUREMENTS AND FINITE-DIFFERENCE NUMERICAL 
CALCULATIONS FOR THREE DIFFERENT SOURCE FIELDS 


Two methods used to study electromagnetic 
induction problems are theoretical model calculations and 
analogue model experiments. Most mathematical methods 
require highly simplified conductivity distributions and 
source fields. Models composed of horizontally layered 
conductors and generalized source fields present no 
serious mathematical difficulties. Problems in which — 
lateral conductivity contrasts occur are much less 
tractable analytically and are usually solved by 
numerical techniques. One such technique which lends 
itself readily to a wide range of three-dimensional 
problems is the finite difference method of section 2.8. 
Various source field configurations can be modelled by 
supplying the appropriate boundary conditions to the 
finite difference mesh. In this manner the electromag- 
netic field solutions for embedded anomalies with the 
Same stratification at all boundaries for various source 
field configurations can be found. Scaled analogue model 
methods are appropriate to problems for which no 


mathematical solution exists or for which numerical 
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evaluation presents great difficulty. This method is most 
readily adapted to problems involving high conductivity 
contrasts (of the order of 107) in which suitably shaped 
graphite bodies simulate highly conducting regions such 
as ore bodies or oceans, and concentrated salt solutions 
Simulate a poorly conducting host earth. A comparison 
of results from the two methods is useful in understanding 
the applicability of both methods as well as gaining an 
appreciation of the difficulties encountered in tnem. 

One difficulty often encountered in laboratory 
scale model experiments is the problem of edge effects. 
Sizeable edge effects can arise in Simulating induction 
problems which involve uniform source fields, since it is 
often difficult to generate perfectly uniform fields in 
the laboratory. Another problem in model experiments is 
that associated with the finite dimensions of the detector 
probes since a finite-sized detector smoothes out the 
behaviour of the electromagnetic fields in regions where 
the fields vary rapidly, i.e.; -in the neighbourhood of a 
Eonductivity discontinuity. There 15 also-o limit to sow 
accurately a finite-sized probe can be positioned. For 
small sized anomalies the positioning and alignment of 
the anomaly is a very critical factor in determining the 
character of anomalous fields. The alignment of the source 


field also becomes more critical as the dimensions of the 
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current carrying device which produces it decreases. 

In numerical modelling techniques a mesh of grid 
points is superimposed on the conductivity configuration 
under study so that the conductivity model is represented 
by a number of cells of uniform conductivity. The 
accuracy of the numerical method depends on the size of 
the grid intervals as well as on the uniformity of grid 
spacings, particularly for points near conductivity 
discontinuities. The volume of the grid array must be 
limited due to storage and computing time considerations. 
This restriction may lead to difficulties in the 
selection of small and uniform grid spacings or in the 
satisfaction of the boundary conditions at great 


distances from the anomaly. 


3.1 Mathematical Analysis 

A brief mathematical development, following 
Dosso (1966), of the modelling problem follows. The MKS 
system of units is used, with the conductivity expressed 
in mho/m. Consider Maxwell's field equations for a 


homogeneous conductor 


3H 
Whe ge A at = 0 (3 1) 
; os ie 
Yael p= ce rie ‘hes 0 (3225) 


To write these equations in a dimensionless form, let 
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where E, H, k._, k_, D, T and S are dimensionless numbers, 
and eo? ho E ial ot df ic ten and O, are the unit quantities 
of the electric field, magnetic field, dielectric constant, 
magnetic permeability, length, time and conductivity, 


respectively. If we substitute (3.3), (3.4) and (3.5) 


waco 43.1) and (3.2), we obtain 
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where v is the dimensionless curl operator taken with 
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Equations (3.6) and (3.7) are dimensionless and have the 
Same form as the original field equations. The solutions 
to these equations are invariant under a change in scale 


if the dimensionless coefficients w, ( and y are invariant. 
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hence, it Tsi.possible! to design a scaled model ‘of. an 
actual geophysical electromagnetic boundary value 
problem. 

If we confine our attention to problems where 
€ and u have the free-space values and let S = 1, 
and eliminate e,/h, from equations (328) to (3210) then 
the necessary and sufficient condition that the solutions 


be invariant under a change of scale is that 


d 
Oe = 
t, = constant or ane constant (sacl 4") 


where Li is frequency and 


d,%, = constant : Saba) 


If we let primed symbols refer to the 
geophysical dimensions and unprimed symbols refer to the 
model dimensions, equations (3.11) and (3.12) can be 


expressed as 
dt =0dettes Wgeeliay) 


do = "dro Con) 
If displacement currents are neglected additional freedom 


im-scaling is obtained. Equation (3.9) can be ignored and 


by combining (3.8) and (3.10) 
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Equation (3.15) indicates the relationship 
that must exist between the conductivities, frequencies 
and lengths involved for the actual problem and its 
scaled model in order that the ratios of the field 
components and phase differences are the same for the 


two problems. 


3.2 Analogue Model Apparatus 

A photograph of the equipment used is shown in 
Plate 3.1. The model consists of an overhead oscillating 
field source (three different ones were used), a large 
tank (244 cm by 168 cm and 76 cm deep) containing 
concentrated salt solution (64 cm deep). of conductivity 
21.0 mho/m simulating the uniform upper layer of a poorly 
conducting earth and a graphite block of conductivity 


pak 10° 


mho/m (see Plate 3.2) which served as an anomaly. 
The graphite anomaly used in the analogue model 

Meas ur eer: and a lucite mount designed and constructed 
by the author to hold it rigidly in the tank are shown 

an. Pratet372.9. To-minimpze the etfects. of the concrete 
floor, the bottom of the tank was lined with a 5 cm thick 
graphite layer of conductivity 1.2 x 10° mho/m. This 

layer of graphite acts as a shield for any remaining 

field at the bottom of the tank. If the tank were much 


deeper (several skin depths), then the effects of the floor 


and earth below would present no problem and the graphite 
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layer would not be necessary. Two stainless-steel sheets 
lining the tank walls (parallel to the beam which 
Supports the measurement probe) and connected by a 
copper wire placed outside the tank permitted the 
currents induced in the salt water to flow parallel to 
the source current (perpendicular to beam). Without 
this arrangement the induced currents for two dimensional 
sources would be constrained to flow perpendicular to 
the source current near the walls of the tank. 

The electric and magnetic field detectors, 
mounted on the ends of lucite rods, were attached to a 
moveable lucite plate held on the rigid beam above the 
tank. The vertical and horizontal magnetic field 
detectors consisted of twin coils 0.10 cm long and 0.64 
cm outer diameter. The coils were designed to nave a low 
inductance in order that the resonant frequency would 
be well above any frequencies of interest. The signal 
from the coil was monitored and amplified by a Tektronix 
502 ,0sctliioscope. The amplitude of the output of the 
oscilloscope was measured using a Hewlett-Packard VITVM 
and the phase angle was measured using a North Atlantic 
Phase Angle Voltmeter. A reference signal was provided 
by a small coil situated at a fixed position near the 
field source. The average field along the surface was 


determined by measuring the voltage difference between 
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points 1.48 cm apart. The electric field detector 
consisted of three probes mounted in a lucite rod with 
the points just protruding from the end and making 
contact with the surface of the salt solution. The two 
outer probes, 1.48 cm apart, were connected by means of a 
two-core shielded cable to the differential input of the 
rextronix 502 oscittoscope. The third probewasituated at 
the mid point between the end probes, was connected to the 
cable shield and hence to ground on the oscilloscope. 

In this way the noise common to both probes was removed. 
The amplitude and phase angle of the amplified signal 
were measured in the same way as described for the 


magnetic field signal. 


3.3 Model Description 

The geophysical model (see Fig. 3.1b), 
considered in this work is a highly conducting body of 
conductivity 0.8 mho/m which is 16 km square and 4 km 
thick embedded at the surface of a poorly conducting 
host earth 126.8 km thick and of conductivity 0.21 x 
10°79 mho/m. The host earth overlays an infinite half 
space of conductivity 1.2 mho/m. Three different source 
field current distributions are used in conjunction with 
the geophysical model; uniform, ye oy and a horizontal 


magnetic dipole. The source frequency in each case is 


OrO Jo Hz. 
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POSITION OF HORIZONTAL DIPOLE AND ye?” SOURCE WITH 
RESPECT TO THE ANOMALY 


E Hty+120)e°"/12) ly(y+120,) 


ANOMALY AND PROFILE POSITIONS 


Figure 3.1 Co-ordinate System, Source and Profile 
Positions 
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The model scaling factors used here are 


ee oy ae 
tia a 
ion a 
Sie tcae 


The primed quantities refer to the natural scale in the 
geophysical problem, while the others refer to the 
analogue model. Using the above conductivity scaling 


9 mho/m 


maetor, graphite of conductivity 0.8 x 10 
represents the highly conducting body, of conductivity of 
0.8 mho/m, while the concentrated salt solution of 

conductivity 21.0 mho/m simulates a poorly conducting host 


: mho/m. The model frequency of 


eh of 0.21 <x 10° 
30 khz simulates a frequency of 0.075 Hz while a scale 
length of 1 cm in the model represents 2 km in the 
natural scale. Thus an anomaly 8 cm square and 2 cm 
thick in the model represents a geophysical anomaly 16 km 
Square and 4 km deep. For this model frequency, the skin 
depths in graphite and in salt water are 0.84 cm and 
63.4 cm, respectively, which correspond to 1.68 km for 
the sea water and 126.8 km for the land. 

The theoretical electromagnetic field solutions 


for the three sources over laterally uniform conductors 


can be calculated for the geophysical model by various 
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methods. For the uniform source the analytical solution, 
of Jones and Price (1970), for a layered earth is used. 
The Fourier series “method of Hibbs and Jones (1976a,b) 

is used to obtain the analytic solution for the yeu 
source. his’ source distribution, where a = 12/720., 

is theoretically approximated by 37 appropriately weighted 
rectangular current distributions of width 10 km and 

these are given in Table 3.1. The rectangular arene 
distribution is approximated by one elemental gaussian of 
half width 10 km as described in Hibbs and Jones (1976d,e). 
The current source flows at a height of 110 km (see 

Fig. 3.la) in the geophysical model and is positioned 
such that the maximum value of the current distribution 
coincides with the center of the anomaly. The analytical 
solution for the horizontal magnetic dipole is calculated 
Using the method of section 2./. The magnetic dipole 
polarized in the negative y direction is placed 17 km 
from the center of the anomaly, see Fig. 3.la, at a 
height of 150 km above the surface. 

The electric field values thus calculated are 
used to provide the boundary and initial conditions for 
the finite difference mesh of section 2.9 and from this 
the electric field associated with the three-dimensional 
embedded conductivity anomaly is determined. The 
anomalous electromagnetic fields due to the conductivity 


structure are assumed to be zero at the boundaries of tne 
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Table 3.1] 


Weighting Coefficients of Rectangular 
Current Distributions 
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finite difference mesh. A 25x25x25 mesh of gr?d' points 


is used in the numerical method. After the electric 
field is determined, the magnetic field quantities are 
obtained from the electric field solution by using 
appropriate finite difference equations. A special 
technique must be used in the theoretical finite- 
difference method when both high conductivity contrast 
e 103) and nigh frequency {>"" 01" Hz)” are us'ed™ in 
conjunction as in this comparison. In the finite 
difference method an arithmetic average of the 
conductivities is normally used and works well in 
obtaining solutions even for high conductivity contrast 
as long as the frequency is not too high. However in 
tne high frequency case a transition zone of conductivity 
must be introduced such that the conductivity contrast 
is never more than 100. In the model this zone of 


mho/m surrounds the anomaly 


EOnductivity 0.21 * 10" 
just inside the anomaly boundary and is one grid spacing 
thick. ‘ This allows the electric field to vary: more 
rapidly across the anomaly boundary, now diffused, and 
More accurately represent the actual field there. 

This technique would also be helpful along the 
Surface boundary wnere the E, caused by the anomaly must 
be taken into account. Since the conductivity is zero 


in the region 2 < 0 an infinite number of transition 


zones would be necessary to allow for a conductivity 
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contrast of no more than 100 at any one boundary. This 
is not practical and was tnerefore not done. However, a 


layer of host earth (o = 0.21 x 107° 


mho/m) one half 
kilometer thick was placed over the conductivity 
configuration to help reduce the effects on ES: This was 
also warranted by tne fact that the analogue model anomaly 
fequires a layer of.salt water (0.1-0.2-cm thick) over 

the anomaly in order to make the electric field measure- 
ments. Since E, is used in the calculation of Hy, an 
effect can be expected on that component which will be 
seen‘in thecsection 3.4. 

In the analogue model the current sources were 
constructed in the following manner. A reasonably 
uniform source field was obtained by using two parallel 
Current carrying wires separated by a distance of 2.4 m 
ana situated at a herghtyof I1v2 m above the surface of 
the salt solution (Ramaswamy et al., 1975). In the model 
Coordinate system, Fig. 3.tbs the current flow in the 
Source=*retd=ts*atong’ thex-direction. The eae source, 
where a = 1./60. (all dimensions in cm), is composed of 
37 current carrying wires, see Fig. 3.2, in which only 
every fifth wire is used. The wires are aligned parallel 
to the x axis and the current carrying elements are 5 cm 
apart. Thus each current element approximates a gaussian 


Current distribution of 10 km half width in the 


geophysical model. 
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The amount of current which flows in any 
particular current carrying element is controlled by a 
resistance placed in series with the wire. The values 
of the resistances vary inversely as the coefficients 
of the rectangular current distributions in the 
geophysical source. The thirteenth wire (the element 
with maximum current) is placed over the center of the 
anomaly. The grid of wires is situated at a height of 
95 cm above the surface of the salt solution. The 
validity of the analogue model for the overhead magnetic 
dipole source has been tested earlier (Thomson et al., 
1972). The magnetic dipole is modelled by a current 
nmoop Of radius 6.4 cm placed at a.height of 75 cm above 
the surface of the salt water solution. This height 
corresponds to a height of 150 km in the geophysical 
model. The current loop was placed at (x = 0 cm, y = -8.5 
cm) with the moment of the dipole polarized in the y 


eryection as shown in Fig. 3.1b. 


3.4 Results 


The parameters which are commonly used to 
describe the effect of conductivity contrasts on the 
geomagnetic field are the ratio of the amplitude of the 
vertical magnetic field (H,) and the component of the 
horizontal magnetic field (H.) that is normal to the 


Source current, and the apparent resistivity (p.) which 
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is defined as 

Pa a lq 

0 y 

calculated in the manner of Cagniard (1953). Profiles of 
these quantities and their relevant phases along with the 
amplitude and phase profiles of the electromagnetic field 
components from which they are derived are calculated in 
the theoretical model and measured in the analogue model 
for 3 positions relative to the anomaly (see Fig. 3.1b) 
for each source field configuration. The profiles of the 
theoretical and analogue models are compared over the 
range -40 km to +40 km in the y direction for the uniform 
source ands ti2 km, for the, other, sources. The, points Tor 
which the calculations are made along the profiles are 
eee Te 0 Zehr 0 5/7 909.5) 0a) 2,16n24.40 Km for, tne. uniform 
Source, and yj= + 0,2.4,6,8,10,12,14,18,26,42 km for the 
other sources. Each field component or ratio with the 
exception of the apparent resistivity and JH 7H | is 
normalized with respect to its value at -40 km or -42 km. 
All phases are normalized to zero at -40 km or -42 km. 
The positions of the profiles, as shown in Fig. 3.1b, are 
across the center of the anomaly, x = O km, four kilometers 
ante the edge of the anomaly, x = -4 km, and outside the 
anomaly four kilometers from the edge for the uniform 
Source, x = -12 km, and ten kilometers from the edge, 


x = -18 km, for the other sources. Also three dimensional 
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figures are given to show the comparison of the electro- 
magnetic field components for the theoretical and analogue 
model. The range of the profiles in the y direction is the 
Same as previously and the positions of the profiles in 
mee x direction are x = 0,-2,-4,+-6,-7,-8,-9,-10,-12,-16 km 
Pom tne uniform source and x = 0,-2,-4,-6,-8,-10,-12,-14, 
=-18,-26 km for the other sources. The figures are dis- 
torted near the edges since a uniform grid spacing is 


used in the plotting. 


3.4.1 Uniform Source Results 

The profiles and three dimensional figures of 
the results for the uniform source comparison are shown 
‘vert os. 3.3 to 3.6. The first set of profiles: given in 
neo 3.3 1S taken 12 km from the center of the anomaly in 
the negative x direction. The normalized electric field 
amplitude, pews of both the analogue model measurements 
and the theoretical calculations exhibits an increase in 
the field due-to the concentration of charge on the 
vertical boundary of the anomaly. The normalized horizon- 
tal magnetic field, JH» of the analogue model is 
reasonably constant across this profile as is the 
theoretical model. The normalized |e curve in the 
theoretical model exhibits a small anomalous field in 
comparison to the analogue model which exhibits almost 


no anomalous field. Since the electric: field in the 
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Figure 3.3 Ge JH I. [HO I> JH O/H | and p, profiles 
with normalized phases for uniform source 
field and profile 1. 
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Figure 3.4 As Figure 3.3 but for profile 2. 
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Figure 3.6 Three-dimensional profiles of [Ele |H\| and 
[HI for uniform source field. 


Theoretical profiles (a,b,c) 


Analogue profiles (d,e,f) 
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theoretical model varies in a less abrupt manner and over 
greater spatial distances than the analogue model electric 
field, it could reasonably be expected that the anoma- 
lous magnetic field of the theoretical model will also 
vary over larger distances than the analogue model 
anomalous magnetic field. This is the case for the [H, | 
field component. Also the magnetic field of the analogue 
model at y = + 40 km is greater than that of the 
theoretical model due to some non-uniformity in the 
analogue model source field, as can be seen in the 

JH S/H | profile. This will tend to lessen any anomalous 
field characteristics in the analogue model curve. The 
JH /HY | profile shows the small amount of non-uniformity 
in the analogue model source field. The apparent 
resistivity curves of both the theoretical and analogue 
models are approximately the same value of 4.76 x 178 
ohm/m, the actual resistivity of the poorly conducting 
host earth. 

The analogue curve exhibits a slight increase 
near the center of the profile, which is due to the 
increased electric field in that region caused by surface 
charge concentration on the anomaly boundary. The phase 
comparisons agree quite well and are nearly constant with 
the exception of the H, and H)/H, phases which shift 


through 180° over the center of the anomaly, which is 
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consistent with the change of sign of a there. 

The second set of profiles, taken along the 
line x = -4 km is given in Fig. 3.4. The effects of the 
anomaly are clearly indicated. rae decreases over the 
anomaly for both the analogue and tneoretical model. The 
theoretical calculation does not vary as rapidly over 
the edges of the anomaly as do the analogue model results 
due to the conductivity transition zone which has been 
introduced at the boundaries of the anomaly of the 
theoretical model. |H | increases over the anomaly for 
the theoretical and analogue models. The theoretical 
model exhibits a greater increase which can be attributed 
to the absence of a transition zone at the surface. 

The inclusion of a transition zone allows the vertical 
electric field E bes more rapidly and thus increase 


the accuracy of the =" term in the calculation of the 


horizontal magnetic field Hy 


|H,| and [H/H | are increased markedly over 
the boundaries of the anomaly for both the theoretical 
and analogue model as would be expected. A spreading of 
the peaks of the theoretical curve is noted for both 
profiles. This is attributable to the slower variation 
of the electric field components in the theoretical model. 
The apparent shifts in the minima of the }H, | and JH 7H | 


Curves are due to the difficulty in positioning the 
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anomalous structure in the analogue model. The comparison 
of the apparent resistivities shows the effect of the 
anomaly. The theoretical and analogue calculations vary 
from the value of the poorly conducting host earth, 
4.76 x 103 ohm/m, to approximately the value of the highly 
conducting anomaly, 0.8 ohm/m. The snape of the theoreti- 
cal curve near the anomaly resembles the shape of the 
pee | curve. The phases of E, and Hy increase over tne 
boundaries of the anomaly whicn can be attributed to the 
transition zone while the phase of P, shows a small 
general increase over the whole anomaly. Again the phases 
of bi and fA, change by 180° as the center of the anomaly 
is crossed. The phase comparisons for this profile 
agree reasonably well. Near the boundaries of the anomaly 
small discrepancies are noted. 

The third set of profiles,-taken across the 
center of the anomaly, x = 0 km, is given in Fig. 3.5. 
The shapes of the amplitude curves in tnis profile are 
much the same as in the second profile. However, the 
PRR) 7 Ps and JH 7H | comparison curves agree more closely 
in this profile. The phase comparisons are similar to 
the last ‘profile with the H, and H)/N) phases slightly 
more accentuated. 

The three dimensional amplitude profiles of the 


comparison of the three electromagnetic field components 


is shown in Fig. 3.6. The bes [H| and JH | three 
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dimensional profiles for the theoretical model are shown 
iieerwds. 3.6a,b,c and for the analogue model in Figs. 3.6 
Geer... Each figure is scaled individually. The increase 
in eee just outside the anomaly due to charge on the 
boundary is clearly evident in the analogue model 
Fig. 3.6d along with the decrease in the field over the 
anomaly. In Fig. 3.6a the smoothing effect of the finite 
Qutrerence calculation can be seen. In Fig. 3.6b and 
3.6e€ an increase in [HI can be seen for both models. 

The erratic behaviour of Fig. 3.6e is due to the small 
variation in the |H | field over the anomaly which when 
combined with the scaling of the three dimensional figure 
Shows the slight errors in manual digitization of the 
analogue model records. The value of |Hy| over the 
boundary of the anomaly parallel to the y axis in Fig. 3.6b 
is higher than would be expected and can be attributed to 
the transition zone. A comparison of Fig. 3.6c and 

Bide. 3.0f shows that [HOI increases over the edges of the 
anomaly for both models. However, the theoretical model 
shows an increase over the boundary of the anomaly parallel 
to the y axis, again due to the transition zone there. 

Tic discrepancies noted in the magnetic field 
components between the analogue and theoretical models 
occur mainly near the transition zone at the edge of the 
anomaly. In the analogue model rapid spatial variations 


due to surface charge distributions are noted in the 
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electric field for profiles just outside the boundaries 
of the anomaly. Such large variations do not occur in 
the theoretical model since the fields vary slowly through 
the transition zone which represents the boundary and the 
charge distributions there are more diffuse and the 
electric fields associated with them are less. Since the 
theoretical electric fields show less effect due to the 
charge in the transition zone, the magnetic field 
solutions are consequently affected and higher values of 
magnetic field are thus generated in the theoretical 
model near the transition zone. The smoothing introduced by 
the transition zone is greater than that of the finite size 
of the sensor in the analogue model since the width of the 


transition zone is twice the probe spacing in the sensor. 
3.4.2 Two-Dimensional Non-Uniform Source Results 

The profiles and three dimensional figures of 
the results for the yer? source comparison are snown in 
mcs 6.7. to. 3. 10.. The first set of profi ies (Fig. 3770 
is taken 18 km from the center of the anomaly in the 
negative x direction. The JE, |» [Hy | and p, comparison 
Curves agree well and are similar in shape to the uniform 
source curves for the profile except for a slight tilting 
due to the non-uniform source field. The |H. | and JH O/H | 
comparison curves show a more definite source effect in 
the form of a general decrease in the amplitude across 
the region of interest. The theoretical and analogue 
models agree well at the y = + 42 km points thereby 
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Fiqunersee heels JH |» [HoI> [HA /H | and p, profiles 


with normalized phases for a two-dimensional 
source field and profile 1. 


- theoretical calculations 


+ analogue model 
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Figure 3.10 Three-dimensional profiles of ee |H.| and 
|H, for the ye °~ source field. 


Theoretical profiles (a,b,c) 


Analogue profiles (d,e,f) 
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Hibbs and Jones (1976d,e). Over the anomaly the 
theoretical curve shows a small anomalous field for 
the same reason as in the uniform source case, while 
the analogue curve does not. The effect of the source 
field is greatest on the [Ho and [HH | profiles for 
this source field configuration. Only a slight source 
effect was noted in the other profiles. The phase 
comparisons agree quite well and are nearly constant for 
Es HY, and Pa: The phases of He and Hie Al vary slightly 
Over the anomaly for the analogue model. 

The second set of profiles, taken along the 
line x = -4 km is given in Fig. 3.8. The effects of 
both the anomaly and the source field are clearly 
indicated... The JE, curves agree quite well with the 
effect of the transition zone clearly evident. The 
source field for both models increases from the negative 
to positive y regions. 

|H| increases over the anomaly and tilts 
upward slightly due to source effect in both the 
theoretical and analogue models. The larger increase in 
|H_ | for the theoretical model occurs for the same 
reason as given earlier for the uniform source. The 
small peaks in the JH | profile at y = +8 km are due to 
the transition zones at those points. |H,| and JH O/H | 


increase over the edges of the anomaly for both 
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theoretical and analogue models though less in the case 
of the analogue model. The peaks are now non-symmetric. 
Both analogue and theoretical curves exhibit the same 
behaviour. The comparison of the apparent resistivities 
Shows the effect of the anomaly but no spatial variation of 
apparent resistivity due to the non-uniform source is evident. 
The phases of E, and p, increase slightly over the boundaries of 
the anomaly for the theoretical calculation. The phase of HY 
remains constant in both analogue and theoretical model. The phase of 
He and Lila dips to -60° at y= -8 km and increases to +60° 
at y = +8 km over the boundaries of the anomaly for the 
analogue model. At y = +42 km the phase of the analogue 
model begins to decrease aS a consequence of the 180° 
phase change which takes place outside the range of 
interest, due to the non-uniform source. The phase 
curves of the theoretical model exhibit approximately 
the same behaviour as the analogue model but due to the 
transition zones the curves do not coincide. 

The third set of profiles taken across the 
gan tem oft the anomaly. x = 0» km, liScagd ven sinpeF iter 3190 
The shapes of the amplitude curves in this profile are 
much the same as in Fig. 3.8 with the exception that the 
JE. | and p, comparison curves agree more closely in this 
profile. The [H\| comparison curves show a greater 


increase over the anomaly than in the previous profile. 
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Both the analogue and theoretical [Ho I> JH O/H, | 
comparison curves are more peaked over the boundaries 
of the anomaly than in the preceeding profile. The 


phases of EY H and Bf Hy are constant for both models. 


y 
The phases of H. and nan are of the same character as 
the previous profiles with an analogue model phase 
decrease to -90° at y = -8 km and increase to +90° 
acy, = +8 km. 

The three dimensional amplitude profiles of 
the comparison of the three electromagnetic field 
components are shown in Fig. 3.10. The |E,|, JH | and |H, | 
three-dimensional profiles for the tneoretical model are 
shown in Fig. 3.10a,b,c and for the analogue model in 
Fig. 3.10d,e,f. Each figure is scaled individually. 
As before the increase in LE | just outside the anomaly 
due to the charge on the boundary is evident in the 
analogue model, Fig. 3.10d, along with the decrease in the 
field over the anomaly. In Fig. 3.10a the smoothing 
erfect of the transition zone can be seen. In Fig. 3.106 
and Fig. 3.10e an increase in |H | can be seen for both 
models. When compared with the three-dimensional figures 
ay 


ene uniform source the effect of the ye source is 


evident in the tilting of the JH | profiles. Tne value 
of [Hy | over the boundary of the anomaly parallel to the 


y axis is higher than would be expected due to the 
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transition zone as in the uniform source case. A 
comparison of Fig. 3.10c and Fig. 3.10f shows that [H, | 
increases over the edges of the anomaly in both profiles. 
The source field causes non-symmetry in the peak values. 


As before an increase in the theoretical les 


figure is 
noted over the boundary of the anomaly parallel to the 


WeraXas.. 


3.4.3 The Horizontal Dipole Source Results 

The profiles and three dimensional figures of 
the results for the horizontal dipole source comparison 
meeesiown. in Fig. 3.11 to Fig..3.74. The first-set: of 
Pearhies, Fig. 3.11, 1s taken 18 km from the center of 
the anomaly in the negative x direction. The concentration 
of charge on the boundary of the anomaly is clearly 
evident in the analogue model electric field JE, | 
Peorileé Also the non-uniform source field and its 
effect on the anomalous electric field is apparent. 
The theoretical curve agrees well with the analogue curve 
confirming the theoretical boundary value calculation of 
Ramaswamy (1973) used here. The magnetic field, ATH) 
profile of the analogue model shows the source field 
Variation but little effect from the anomaly. The 
theoretical field curve increases from negative to 
positive y but by an insufficient amount to exactly match 


the analogue curve. This misalignment of the curve is 
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Figure. 0 | lbeai JH | eae JH O/H and 9, profiles 
with normalized phases for the horizontal 
dipole source field and profile l. 
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Figure 3.13 As Figure 3.11 but for profile 3. 
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Figure 3.14 Three-dimensional profiles of JE, I> }H | anc 
|H, | for the horizontal dipole source field. 


Theoretical profiles (a,b,c) 


Analogue profiles (d,e,f) 
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A? 


WH, 
explained by the fact that the value of the magnetic field 
at y= -42 km used to normalize the curve is Slightly less 
accurate due to the wider grid Spacing used at that point. 
This, combined with the larger spatial variation of the 
source field, serves to accentuate the misalignment when 
compared with the [Hy | profiles of other sources. The hale 
and [HH | profiles compare quite well with this profile. 
The theoretical curves show a slightly more anomalous field 
as was the case for the uniform source. The apparent resis- 
tivity has the value of the poorly conducting host earth, 


4.76x 10° 


ohm/m, for both models. The analogue model curve 
shows a slight increase over the anomaly due to the 
increased electric field. This increase is not as large as 
noted for previous sources for this profile. This is a 
source effect due to the fact that [Hy | has a larger value 
over the anomaly for the horizontal dipole source. The 
theoretical curve is slightly higher than the analogue 
Curve due to the slightly lower value of the theoretical 
[Hy | curve. The phase comparisons agree reasonably well 
and are nearly constant with the exception of the H, and 


H_/H. phases which shift through 180° at the minima of 
at 

the amplitude curves of these quantities. | 
The second set of profiles taken along the line 
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of the anomaly are clearly evident in this profile. Both 


the analogue and theoretical JE,| curves decrease over the 
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V7 
anomaly by approximately the same amount and exhibit the 
Same source field with the exception of the effect of the 
transition zone. 

Near y = +42 km both the theoretical and 
analogue |H | curves assume the same relationship as 
in the previous profile, however, over the anomaly both 
curves now increase. The theoretical curve tilts 
upward in the positive y direction over the anomaly with 
the two small peaks at y = +8 km. The analogue curve 
exhibits a decrease over the edge of the anomaly at 
y = -8 km and an increase over the edge at y = +8 km. 

The |H,| and JH O/H | curves for the theoretical 
and analogue models exhibit two non-symmetric peaks over 
the edges of the anomaly along with the profile of the 
source Raita The theoretical curve peaks are slightly 
higher and wider as was the case for the uniform source. 
The theoretical and analogue apparent resistivity curves 


vary from approximately the poorly conducting value of the 
3 


host earth, 4.76 «x 10> ohm/m, ‘to the value of the highly 


conducting anomaly, 0.8 ohm/m. The phases of EY and Pa 
increase slightly over the boundaries of the anomaly for 
the theoretical calculation which can be attributed to the 
transition zone at the boundaries. The phase of HY 
remains constant in both the analogue and theoretical 


models. The phase of H, and HSH dips to -60° at 
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y = -8 km and increases to +60° at y = +8 km over the 
boundaries of the anomaly for the analogue model. A 180° 
phase shift is noted at y = +18 km due to the source 
field in the analogue model. The phase curves of the 
theoretical model exhibit first an increase and then a 
decrease at y = +8 km due to the transition zone. The 
theoretical curve then remains constant until the 180° 
phase shift at y = +18 km. 

The third set of profiles taken across the 
center of the anomaly, x = 0 km, is given in Fig. 2.18. 
As was the case for the previous sources these profiles 
have much the same shapes for the amplitude curves as 
the previous profiles of Fig. 3.12. The een and p, 
comparison curves more closely agree in these profiles. 


|. |H 


The trends of the |H and JH /H | curves seen in 


y z| 


the previous profiles for the analogue model are fully 
Berined in this profile. The increase of the. theoretical 
|H | curve is more pronounced. The phase profiles are 
much the same as in the previous profile. 

The three dimensional amplitude profiles of the 
comparison of the three electromagnetic field components 
fee chown in Fig. 3.14. The Rese |H | and |H. | three- 
dimensional profiles for the theoretical model are shown 
in) Fig. 34.l4a,b,c and for the analogue model in Fig. 3.14 


Gee.f and-each figure is scaled individually as before. 
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As with the other sources the increase in [EY | just 
outside the anomaly due to the charge on the boundary is 
evident in the analogue model, Fig. 3.14d, along with 
the decrease in the field over the anomaly. In Fig. 
3.14a the smoothing effect of the transition zone can be 
seen. The source effect is clearly evident in the JH | 
and |H, | profiles for the analogue model, Figs. 3.14c,f. 
The values of the theoretical |H | and [Ho | curves, 
Fig. 3.14b,c, over the boundary of the anomaly parallel 
to the y axis are higher than would be expected due to 
the transition zone. A comparison of Fig. 3.14c and 
Fig. 3.14f shows that |H. | increases over the edges of 
the anomaly for both the theoretical model and analogue 
model. However, the analogue model variation of |H, | is 


less accentuated. 


3.5 Summary 

The two very different methods for studying 
electromagnetic perturbation problems show general 
agreement. The condition of a gradually changing 
electric field through transition zones in the theoretical 
calculations leads to differences in the results near 
the boundary of the anomaly especially for the magnetic 


field. This effect, however, has been lessened by the 


inclusion of the transition zone. This effect Vs both a 


function of the frequency, the conductivity contrast and 
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to a lesser extent the grid spacing. The solution over 
the anomaly and at large distances agrees well with the 
analogue model thereby confirming the boundary condition 
calculations of Hibbs and Jones (1976d,e) and 
Ramaswamy (1975). The frequency, .075 Hz, together 
with the high conductivity contrast tests the limits of 
the approximations in the methods. It is seen that a 
transition zone must be introduced in the numerical 
method to allow the fields to vary more rapidly. For 
longer periods or lower conductivity contrasts the 
Biel ds vary less rapidly and are better approximated 
particularly in the numerical calculations. It has thus 
been instructive to examine the possible limits on the 
approximations applied to the iterative process due to the 
complexity of the process and the large number of 


configurations which can be modelled. 
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CHAPTER 4 


A STUDY OF PERTURBATION AND INDUCTION ARROWS FOR 
VARIOUS CONDUCTIVITY CONFIGURATIONS AND SOURCE FIELDS 


The numerical methods of Chapter 2 are used 
to calculate the electromagnetic fields associated with 
three dimensional conductivity anomalies. The calculations 
were done for three source fields and three conductivity 
configurations. A uniform field oriented in various 
directions relative to the anomaly is considered first. 

The transfer functions and related perturbation and 
induction arrows associated with the fields are calculated. 
The results show that the perturbation arrows do not well 
describe tne anomalous currents for the three-dimensional 
model considered, but can give a method of outlining tne 
spatial extent of the anomaly. The induction arrows may 

be used in the traditional way to indicate the flow of 
| anomalous currents and point toward conductivity 
inhomogeneities. Two methods of calculating the induction 
arrows are compared. 

The second type of source field considered is 
due to a number of two-dimensional sheet current sources 
with different current intensity distributions, variously 
positioned and orientated witn respect to the conductivity 


anomaly. The transfer functions and related perturbation 
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and induction arrows associated with the fields are 
calculated as before. In addition, approximate normal and 
perturbation fields are used in the transfer function 
calculations and comparisons are made. 

The third source field considered is due to 
horizontal magnetic dipoles placed at two different 
positions with respect to the conductivity anomaly. The 
transfer functions and related perturbation and induction 
arrows associated with the fields are calculated and 
compared with the arrows obtained from a uniform source 
calculation. Again, the results show the source effect on 
the induction arrows and indicate that the perturbation 
arrows provide a method of outlining the spatial extent of 
the anomaly. The transfer function calculations are made 
for both exact and approximate normal fields. In the 
transfer function calculation the anomalous fields are 
correlated with a normal field as suggested by Schmucker 


(1970) and Cochrane and Hyndman (1970). 


4.1 The Mathematical Formulation of the Induction and 
Perturbation Arrows 
In studies of geomagnetic variations, it is 
common to seek empirical relationships between different 
magnetic field components which are independent of time. 


The electromagnetic induction phenomena encountered in 


geomagnetic studies obey Maxwell's equations. MTnese 
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equations are linear, and thus linear induction matrices 
can be defined. The determination of the elements of these 
matrices and from them the conductivity structure of the 
earth is the goal of transfer function analysis. Schmucker 
(1970) defined two sets of arrows, the induction and 
perturbation arrows, in terms of the elements of his 
maauctiron matrix. 

In his approach, Schmucker (1970) assumed that a 
linear relationship exists between the Fourier transforms 
of the anomalous internal field components and the normal 
field components for a single frequency and this is 


expressed as: 


C h h h E § 
x y Zz Xn Ax 
= d c oa C oO 
"Wa dy dy Zz Vi Ay ( ) 
C 4 z Z C C 
Z 5 X y Z Ze AZ 
iti and  € are Fourier 
The quantities oe : a n a 


transforms of the anomalous parts of the x, y and z 


components of the magnetic field respectively; a , Cy 


ang C are the Fourier transforms of the normal parts 
z 


n . 
of the x,y and z magnetic field components respectively; 


n 


are the Fourier transforms of the 
and C, 5 Cry and C,, 


uncorrelated portions of the x, y and z components, where 
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the Fourier transform for a function R(t) .iS pdefinedas 


+00 
SACP) ss c(phGt ie wer oe (4.2) 
puewtransfer funttions h, d and.z.aré in general complex 
and their real parts are referred to as the in-phase 
transfer functions while their imaginary parts are the 
quadrature-phase transfer functions. Then for the 


anomalous part of Z. 


G yy ee @ er Sala & sy eae + C (4.3) 


The autocorrelation power spectrum of C5 is defined as 


a 
Cc oma (4.4) 
S = ° , 
Lee Z. Z. ) 
where ie is the period chosen, Sz z ts the Vinpting 
aa 


value as Te > o and oF denotes the complex conjugate of vom 


In a similar manner the cross correlation power spectrum 


between C5 and Cy for example is defined as 
a n 


(445) 


Using the crosscorrelation and autocorrelation 
power spectra of the anomalous and normal field components 
along with the condition that the crosscorrelation power 


Spectra between the residual and normal field components 
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be zero, Schmucker (1970) obtained a linear matrix 
involving the transfer functions and power spectra. For 


the z transfer functions: 


Z oe 7d Ly es aS 
Xn Y ¥nXn 4 7%, ZX, 
205 +03RS Feet S = § (4.6) 
X Xan y YnXn Z 2nYn Zan 
zUs + z § +275 = § 
Xan y Yn7n a zn 7 n Zacn 


Similar matrices can be obtained for the h and d transfer 
functions. 

Schmucker (1970) represented the transfer 
functions as two sets of arrows. If we let i and j be 
Cartesian unit vectors in the x and y directions respective- 
ly, the induction arrows are defined for a particular 


frequency as: 


C eae i Sez j (in-phase arrow) 
~real Xreal Yreal 

aD) 
ae . a7 qeeteez. j= (quadrature- 
imaginary Ximaginary — Yimaginary phase arrow) 


in accordance with the orientation of the Parkinson (1959) 
arrow. The in-phase arrows, as defined here, point toward 


Current concentrations and thus toward regions of higher 


conductivity than the surroundings. 
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The second set of arrows, called perturbation 


arrows, is defined as: 


Ddpelyhet ydyd 
(4.8) 

qi pet dd 
for both in-phase and quadrature-phase arrows. Schmucker 
(1970) points out that when these arrows are rotated 90° 
counter-clockwise they indicate the strength and direction 
of the anomalous internal current field which is super- 
imposed on the unperturbed normal current flow. The p 
arrow represents the anomalous current in the - y 
direction and the gq arrow that in the x direction. 

In the present work the in-phase induction 
arrows are made negative with respect to the definition 
by Schmucker (1970) and consequently point away from high 
conductivity areas toward areas of low conductivity. The 
perturbation arrows are rotated 90° counter-clockwise in 


the manner suggested by Schmucker (1970). 


4.2 The Calculation of Perturbation and Induction Arrows 
for a Uniform Source Over a Three-Dimensional 
Conductivity Model 


The calculations are carried out for uniform 


Sources linearly polarized at various directions over a 


conductivity inhomogeneity embedded in a two-layered semi- 


infinite conducting region with a plane boundary. The 
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first layer is 4 km in depth and of conductivity 4 mho/m. 
The second layer is of infinite depth and has conductivity 
0.001 mho/m. The embedded anomaly is a 6 km square intru- 
sion of the lower conductive layer which penetrates the 
top layer to the surface of the conductor. 

This conductivity configuration was chosen 
because of its symmetry since all the electric field 
solutions for the source polarization angles presented 
could be derived from three independent solutions (0°, 
Fe aS? ).. 

The numerical technique of Section 2.8 is used 
to obtain the electric field solution from which the 
magnetic fields and transfer functions are calculated. 

The calculations have been made for two frequencies: 


Pev0) Hz and 0.0033 Hz. 


4.2.1 Single Source Results 

For linear polarization and a single orientation 
of the source field in which the normal magnetic field is 
in the y direction the transfer function matrix reduces 


to a single equation. For the z transfer function matrix 


this equation is: 


Similar equations can be obtained from the h and d 


transfer function matrices. 
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The transfer functions obtained in this manner 
for the model considered are shown in the arrow representa- 
tion in Figure 4.1 for 0.001 Hz. Since the normal current 
is polarized in the x direction, the p perturbation arrow 
is undefined. The gq perturbation arrow should then describe 
the magnitude and direction of the anomalous current. The 
q arrow together with the arrows associated with the 
calculated anomalous current (from the known electric 
perturbation field and the conductivity structure) and the 
induction arrows are shown in Figure 4.1. 

Away from the anomaly the perturbation arrows 
(Figure 4.1A) qualitatively resemble the anomalous current 
arrows (Figure 4.1B). Nearer the anomaly the in-phase 
perturbation arrows differ markedly from the in-phase 
anomalous current arrows. Over the anomaly the in-phase 
perturbation arrows greatly increase in magnitude. Since 
the conductivity of the anomaly decreases by a factor of 
4 x 10°, the current flowing in the anomaly is small 
relative to that flowing in the surrounding material and 
therefore only small anomalous current arrows appear over 
the anomaly. In the method of Section 2.8 the conducti- 
vity associated with a particular node is the average of 
the conductivities of the surrounding regions, and the 
boundary between the anomaly and the Surrounding region 
is not abrupt as indicated in the figure but consists of 


a transition zone between the conductive regions. 
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Figure 4.1 Single source with the normal current field in 
the X direction for 0.001 Hz. The dashed 
arrows are the in-phase arrows and the solid 
arrows are the quadrature-phase arrows. 


(A) q perturbation arrows, 
(B) anomalous current arrows, 
(C) induction arrows. 


The surface of the conducting region near the 
anomaly is shown and the inner square region 
represents the anomalous structure. 
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Figure 4.1C gives the induction arrows associated 
with this simple source orientation. The arrows are 
polarized in the direction of the normal magnetic field. 
They vary in magnitude over the grid and point toward the 


conductivity anomaly. 


peewee) Multiple Source Results 

For a number of sources with the inducing field 
linearly polarized and oriented in various different 
directions the transfer function matrix is represented by 


a two-dimensional matrix. For the z transfer functions 


(4.10) 


The cross correlation and autocorrelation power spectra 
are calculated for each source polarization and these are 
averaged and the transfer functions are determined in the 
manner of Schmucker (1970). 

The results from the calculation of transfer 
functions for six sources for periods of 0.001 Hz and 
0.0033 Hz are shown in Figures 4.2 and 4.3. For the case 


of multiple sources at different orientations both p and 


q arrows are defined. By adding these two arrows vector- 


ially (p + q) an arrow to represent the total anomalous 
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Figure 4.2 


Perturbation, current and induction arrows 
calculated from six source orientations: 23°, 
45°, 90°, 135°, 248° and 315° as measured 
counter-clockwise from the x-axis for 0.001 az: 
The dashed arrows are the in-phase arrows 

and the solid arrows are the quadrature-phase 
arrows. 


(A) p perturbation arrows ; 

(B) q perturbation arrows ;» 

(C) anomalous current arrows ;» 
(D) p + gq perturbation arrows ; 
(E) induction arrows. 
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Figure 4.3 As Figure 2 buts for (0.00338 .nzZ- 
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current should be obtained. The p and q arrows together 
with the arrows representing p + q are given in these 
figures. Also, vectors which describe the calculated 
anomalous electric currents determined from the anomalous 
electric field and known conductivities are shown in 
Figures 4.2C and 4.3C. It is apparent that the (p + q) 
arrows do not compare well with the anomalous current 
arrows for this three-dimensional model. From Figures 
4.2— and 4.3E, we see that the induction arrows are no 
longer polarized, but now point toward the center of the 
anomaly over the whole surface region. The stability of 
the perturbation and induction arrows with frequency for 
this model is apparent by comparison of Figures 4.2 and 4.3. 
If the p perturbation arrow is made negative 
with respect to the definition by Schmucker (that is, it 
now defines the anomalous internal current superimposed on 
the normal currents flowing in the + y direction) as in 
Figure 4.4A and then added to the q perturbation arrow 
to give a new (p + q) arrow as in Figure 4.4E, this new 
(p + q) arrow now agrees qualitatively with the calculated 
anomalous current arrows (Figure 4.4D) for regions away 
from the anomaly. As the anomaly is approached both the 
in-phase and quadrature-phase (p + q) arrows differ from 
the anomalous current arrows. However, both sets of 
arrows exhibit the same general trend. Over the edge of 
the anomaly the (p + q) perturbation arrows point approxi- 


mately opposite in direction to the anomalous current arrows. 
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Figure 4.4 


Perturbation, electric field, current and 
induction arrows calculated from six source 
orientations: 23°; 45°, 90°%. 135”, 248° and 
315° at 0.001 Hz. The dashed arrows are the 
in-phase arrows and the solid arrows are the 
quadrature-phase arrows. 


(A) p perturbation arrows, 

(B) q perturbation arrows, 

(C) anomalous electric field arrows, 
(D) anomalous current arrows, 

(E) p + q perturbation arrows, 

(F) induction arrows. 
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Schmucker (1970) stated that at least five random 
source field orientations are required to assure stable 
transfer functions. Figure 4.5 gives a comparison of 
induction arrows for three sets of polarization angles at 
0.001 Hz. In Figure 4.5A two source orientations which are 
mutually perpendicular are used and the induction arrows 
compare well with those of Weidelt (1975) calculated by 
Superposition for a somewhat similar model. In Figure 
4.5B three source field orientations spread over a range 
of 45° are used to construct the arrows and some skewing 
of the induction arrows is evident. Figure 4.5C is for 
six source field orientations which are spread throughout 
tne full 360° range and in this there is very little 
skewing of the arrows. It appears that less than six 
source orientations may be used, but if this is done the 
direction of the source fields must be chosen with care 
and should be quite different from one another. 

Cochrane and Hyndman (1970) proposed that the 
total magnetic field be substituted for the normal magnetic 
field in the transfer function calculations since this 
would be easier to apply in the analysis of experimental 
observations. Figure 4.6 is a recalculation of the arrows 
as in Figure 4.2 but with this change, and it is seen that 
no appreciable change in the configuration of the arrows 


is produced. This is expected because the horizontal 
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Figure 4.5 


Comparison of 


number and angular spread of 


current sources. The dashed arrows are the 
in-phase arrows and the solid arrows are the 
quadrature-phase arrows. 


(A) induction 
angles of 
(B) induction 
angles of 
(C) induction 


angles of 
CU eonne 


arrows for source polarization 
0° and 90°, 


arrows for source polarization 
0°. 23° “and 45 .., 


arrows for source polarization 
23°, 45°, 90°, 135°. 2465 Bane 
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Figure 4.6 


Perturbation, current and induction arrows 
calculated as in Figure 4.2, but with the 
substitution as suggested by Cochrane and 
Hyndman (1970)... AS Bs Ca aD and ue pages 

as. in Pigurne 4.2: 
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magnetic field decreases by only about 25 percent over 
the anomaly for this model. 

In Figures 4.4C and 4.4E a comparison of the 
anomalous electric field and the new (p + q) perturbation 
arrows is made for six source orientations at 0.001 Hz. 
Both the anomalous electric field and perturbation arrows 
are sensitive to the anomaly. The electric field is more 
sensitive to anomalous structure than the magnetic field. 
It may therefore be desirable, since in most experimental 
work of this nature the magnetic fields are measured, to 
obtain electric field information from quantities 
calculated from the observed magnetic fields. In this 
respect, the (p + q) perturbation arrow should be useful 


in accurately outlining regions of anomalous conductivity. 


4.3 The Calculation of Perturbation and Induction Arrows 


for a Three-Dimensional Conductivity Structure and 


Various Two-Dimensional Source Fields 

The calculations were carried out for eight two- 
dimensional sheet current sources with spatial intensity 
as in Figure 4.7 oscillating with a circular frequency of 
2nt/w and flowing at a height of 110 km above a two-layered 
semi-infinite conducting region with a plane boundary and 
an embedded conductivity anomaly. The first layer is 50 
km in depth and of conductivity 0.01 mho/m. The second 


layer extends to infinity and has a conductivity of 0.1 
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mho/m. The embedded anomaly is L-shaped, as in Figure 4.7 
and has a depth of 3 km. The conductivity of the anomaly 
is 0.0001 mho/m, and the calculation is done for a 
frequency of 0.1 Hz. 

The boundary conditions for the finite difference 
mesh were determined by using the general method for two 
dimensional sheet current distributions over layered 
conductors as described in Section 2.6. Each source was 
approximated in a piecewise linear manner by 49 suitably 
weighted rectangular current distributions of 30 km width. 
The rectangular current distributions were approximated by 
Pree elemental gaussians of 10 km half-width each 
Spatially shifted one half-width. These field values 
were then substituted into the numerical technique of 
Section 2.8 and the electromagnetic field solutions 


obtained. 


4.371 Single Source Results 

The normal magnetic field for a non-uniform sheet 
current source polarized in the x direction will contain 
magnetic field components in the y and z direction. 
However, the z component is small in comparison with the 
y component over the anomaly and surrounding region. For 
this reason the 3x3 transfer function matrix (equation 4.6) 
is ill-conditioned and must be reduced to a one-dimensional 


relationship. For the z transfer function matrix this 
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equation is: 


zs CAST) 


Similar equations are obtained from the h and d 
transfer function matrices upon reduction to the one- 
dimensional case. 

The non-uniform sheet current source with the 
electric field polarized in the x direction used in this 
model is shown in Figure 4.7A. The transfer functions 
obtained in this manner for the model considered are shown 
in the arrow representation in Figure 4.8. Since only hy 


ds; and z, are defined for this source, the p perturbation 


y? 
arrow is undefined. The gq perturbation arrows (Figure 
4.8A) should then describe the magnitude and direction of 
the anomalous current. At some distance from the anomaly 
as was found in Section 4.2, the q perturbation arrows 
qualitatively resemble the anomalous current arrows. 
However, over the anomaly the perturbation arrows increase 
in magnitude and are not representative of the anomalous 
Currents. This is to be expected since the anomalous 
Magnetic field is most intense over the anomaly. Also, 
for this source, the induction arrows, Figure 4.8B, are 
polarized in the y direction. The induction arrows vary 
in magnitude over the grid and point toward the anomaly. 


As would be expected the null of the induction arrows 
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Figure 4.8 


Single source with normal current field in 

the x direction for 0.1] Hz. The dashed arrows 
are the in-phase arrows, the solid arrows are 
the quadrature-phase arrows. 


(A) gq perturbation arrows, 
(B) induction arrows. 
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4.3.2 Multiple Source Results 

For a number of sources with inducing fields - 
linearly polarized and orientated in various directions 
the transfer function matrix must now be represented by a 


two-dimensional matrix. For the z transfer functions: 


(4.12) 


Similar equations are obtained from the h and d transfer 
Fiction matrices. The cross correlation and auto- 
correlation power spectra are calculated for each source 
polarization and these are averaged and the transfer 
functions are determined. 

The two sources used to represent a multiple 
source orientation are shown in Figure 4.7B,C. An approximate 
separation of 70° in the orientation of the source 
polarizations allows the calculation of the transfer 
functions with less than six source field polarizations. 

The results of this calculation for the frequency 


used are shown in Figure 4.9. Since h and d transfer 
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Figure 4.9 


Two sources with polarizations separated by 70° 
situated near anomaly for 0.1 Hz. The dashed 
arrows are the in-phase arrows and the solid 
arrows are the quadrature-phase arrows. 


(a) p perturbation arrows, 

(b) g perturbation arrows, 

(c) p + g perturbation arrows, 
(d) induction arrows. 
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25 
function matrices are two-dimensional in x and y both the 


p (Figure 4.9a) and q (Figure 4.9b) arrows are defined. 
The p perturbation arrows have been made negative with 
respect to the definition of Schmucker (1970) (that is, 
they now define the anomalous internal current superimposed 
on the normal currents flowing in the + y direction. If 
these two arrows are added vectorially (p + q) as in 
Figure 4.9c, the spatial extent of the anomaly is outlined 
by arrows polarized at an angle of approximately 135° 
counter-clockwise to the positive x axis. The (p + q) 
arrows are approximately 50% larger than the surrounding 
arrows over the anomaly. 

In this case the induction arrows (Figure 4.9d) 
are no longer polarized in the y direction. The arrows 
tend to point radially toward the central moment of the 
anomaly. However, some source effect due to the non- 
uniformity of the source field can be noted in the in-phase 
arrows. 

As pointed out by Cochrane and Hyndman (Poy Oh 

the problem of defining the normal field for the calcula- 
tion of induction arrows can be partially overcome by 
substituting the total field as the normal field in the 
transfer function calculations, since this is easier to 
apply in the analysis of experimental observations. In 
their calculations Cochrane and Hyndman (1970) assumed 


that the magnetic field component in the z direction was 
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216 
totally anomalous which thus allows a calculation of the 
induction arrows. This assumption cannot be made when 
calculating the anomalous field for the horizontal compo- 
nents of the magnetic field. Figure 4.10 is a recalcula- 
tion of the arrows in Figure 4.9 using the Cochrane and 
_ Hyndman (1970) suggestion. As can be seen in Figure 4.10, 
little change is noted in the perturbation and induction 
arrows, which depend on the horizontal magnetic field 
components. However, the anomalous field used to calculate 
the perturbation and induction arrows in Figure 4.10 is 
exact. 

It is important to determine the effect on the 
calculations if the anomalous and normal fields are 
approximate. In Figure 4.11 the perturbation and induction 
arrows are given for a calculation in which the normal 
field is assumed to be constant and the anomalous field is 
the result of subtracting the normal field from the total 
Magnetic field. The value for the normal field was 
obtained by averaging the field values along the boundary 
of the finite difference mesh. The boundary forms a square 


with a perpendicular distance of 58 km from the center 


of the anomaly. Figure 4.11 gives an indication of the 


relative stability of the Schmucker method. 


Schmucker (1970) stated that at least five ran- 


dom source field orientations are required to assure 
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Figure 4.10 As Figure 4.9 but with substitution as 
suggested by Cochrane and Hyndman (1970). 
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Figure 4.11 As Figure 4.9 but for averaged normal field. 
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Stable transfer functions. Figure 4.12 shows the result 
of calculating the induction arrows for six sources 
centered a large distance from the anomaly with an angular 
spread of 63°. The sources used in this calculation are 
shown in Figure 4.7a, d, e, f, 9, h. In Figure 4.12a the 
calculation was done in the manner of Schmucker (1970) and 
it is apparent that severe skewing of the induction arrows 
is evident when compared to Figure 4.9d. The calculation 
of the induction arrows according to Cochrane and Hyndman 
(1970) is shown in Figure 4.12. Again severe skewing of 
the arrows is noted. However, the appearance of the two 
diagrams differs much more than was the case between 


Figure 4.9d and Figure 4.10d. 


4.4 The Calculation of Perturbation and Induction Arrows 


for a Three-Dimensional Conductivity Model and Dipole 


source Fields 
The calculations were carried out for horizontal 
magnetic dipole sources indicated by |] and 2 as in Figure 
4.13 oscillating with a circular frequency of 27/w and 
placed at a height of 150 km above a two-layered semi- 
infinite conducting region with a plane boundary and an 


embedded conductivity anomaly. The first layer is 128 km 


-3 
in depth and of conductivity 0.21 x 10 mho/m. The 


second layer is of infinite depth and has conductivity 


0.8 mho/m. The embedded anomaly (indicated by the inner 


' 
i 
a 
iy be 


J “1 p=) Ue a ok Dirt es a Bae 
ai oe ; he ad eo ae 
tie 

os ® 

Fy 
ice ats 


_ oe) 
> 
vr. 
= ae 
hy 


tiueesy ott ewonde ST... o1up ta, -anotsonut vorencd 
earivoee xt2e vo? 2worts m0 rs ubAT sits entsstusTa 


at 
VYefupns o6 Adtw Viameas ont mont annede rb’ aptet 6 be 


a16 nofrisivolso~2ras. we baey 2921002 ont | Ty ¥o" se 
oft oST.8 syuprt ol «8 8° AP ae gb BM. b sue ata 

bag (OVOT) vodtoumeAse +0 1oAnbm SAD at wai 25W nobis 9 
eworrs MOTFIUHWT sad to parwote 19" 9% teas norsnas & 7 
aottal ate on? . 68.8 supa ad bavsqmo> new dnab ty 
némbayh bs sfsidood oF on Pb10d26 | aWwovts nottoubnt 
to priwote ovsvee atapA .SF,8 owet ot awore at ( 
owt-9d¢ Yo sonsysoqas ond -¥svewoH ..bston 2t a 
ngewied seso sat esw-neds siom doum) erat tb. ems" 


bors hid si be. ¢ 


[etwostyed vot suo Ce ee 
omup i ar 25 &§ bos T yd beteatbnt 2991002 otogts 3 
bite wNag. do yonewpey? walworta 6 dtiw anisett 20. 
-imee baveysl-ows 6 svods mA OG! to tdpted 6 “tet 

* 7B bine wisbauad snelq 6 site: notes enttoubnos osha 
ma asi et +9MET tert? oft vlemons wtvisoubnes be 
eiT .m\odm “~ Or x FS. 0. vitvtsoubmoo to bis wan 0 
\yJivitsybaos 26f bas dtqeb ot inttat to at newt | 
yonnt oft yd bodsotoat) ulstons bebbodms on 


baeterh Kre 4, ewer) no hess ob ia, 9% 
788) bedsqe slupiis bedla f . | 


Wi eles umn Se ot dafbyoo 33% ‘hotseraates-t > a 
Benen Ot Po fb1b2568 noftstTusts ta) “ Le 


, a a - —") inal » = oa De ee Se ee _ 


\ 


ave . ~ be ome pit _ - 
tT; r. fs My Ps yy ye : ' 
Vo A arr 


_ 
‘ PP AA UF oe 


Figure 4.12 


Induction arrows for six distant sources 
with close angular spread (68°). 


(a) calculation according to Schmucker (1970), 
(b) calculation according to Hyndman (1970). 
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Figure 4.13 Coordinate system, anomaly and two dipole 
sources (not to scale). 
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square of Figure 4.13) is 16 km square and 4 km deep with 

its top at the surface of the conducting region and has 

the same conductivity as the lower layer. With the Origin 

of the coordinate system at the center of the anomaly, 

the position (in Cartesian coordinates) of the first dipole 

is (x=0 km, y=-17 km, z=-150 km) and of the second is 

(x=-17 km, y=0 km, z=-150 km). The dipole moment of the 

first dipole is polarized in the negative y direction and 

the second dipole is polarized in the negative x direction. 
The numerical technique of Section 2.8 

is used to obtain the electric field solution for the 

three-dimensional model and from this the magnetic fields 

and transfer functions are calculated. The boundary 

conditions for the finite difference mesh are obtained 

for the dipole source above the layered conductor by 

means of a program supplied by V. Ramaswamy, (Ramaswamy, 

1973). Near the boundary the anomalous fields are 


asSumed-to be zero. The calculations were made for the 


frequency 0.075 Hz. 


4.4.1 Single Source Results 


For a single magnetic dipole at position 1 the 
theoretically calculated normal magnetic fields contain 


all three field components. However, the x and z 
components are small in comparison to the y component 


Over the anomaly and surrounding region due to the 
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position of the dipole at x=0 km, y=-17 km and z=-150 km. 
For this reason the 3x3 transfer function (equation 4.6) 

is ill-conditioned and must be reduced to a one-dimensional 
relationship. For the z transfer function matrix this 


equation is: 


Zarate (4.13) 


Similar equations are obtained from the h and d transfer 
function matrices upon reduction to the one-dimensional 
case. 

The transfer functions obtained in this manner 
for the model considered are shown in the arrow representa- 


Pon in Figure 4.14. ‘Since only h,., d. and Zy are defined 


yey. 

for this source, the p perturbation arrow is undefined. 
The q perturbation arrows (Figure 4.14A) should then 
describe the magnitude and direction of the anomalous 
Current. At some distance from the anomaly as was found 
previously, the g perturbation arrows qualitatively 
resemble the anomalous current arrows. However, over the 
anomaly the perturbation arrows increase in magnitude and 
are not representative of the anomalous currents. Also, 
for this source, the induction arrows, Figure 4.14B, are 


polarized in the y direction. The induction arrows vary 


in magnitude over the grid and point away from the anomaly. 
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Figure 4.14 Single horizontal dipole source with moment 
polarized in y direction positioned at 
x=0 km, y=17 km and z=150 km. The dashed 
arrows are the in-phase arrows and the 
solid arrows are the quadrature-phase arrows. 


(A) g perturbation arrows, 
(B) induction arrows. 


The surface of the conducting region near 
the anomaly is shown and the inner square 
region represents the anomalous structure. 
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For purposes of comparison the q perturbation 
and induction arrows for a uniform source with a normal 
magnetic field polarized in the + y direction are given 
in Figure 4.15. In the g perturbation arrows the 
quadrature-phase arrows are greater in magnitude compared 
with the in-phase arrows for the uniform field source. 
Also, the q perturbation arrows are more perturbed in the 


uniform field case. The induction arrows are very similar. 


4.4.2 Two-Dipole Source Results 

For the two-dipole source the magnetic field 
solution of dipole 1 was rotated 90° to obtain the 
solution for dipole 2. The cross correlation and auto- 
correlation power spectra are calculated for each indivi- 
dual dipole source and these are averaged and the transfer 
functions determined. The transfer function matrix must 


now be represented by a two-dimensional matrix. For the 


Z aranster functions: 


(4.14) 


Similar equations are obtained from the h and d transfer 


function matrices. 
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Figure 4.15 As Figure 4.14 but for a uniform source with 
a normal magnetic field polarized in the y 
direction. 
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The results of this calculation for the frequency 
used are shown in Figure 4.16. Since h and d transfer 
function matrices are two-dimensional in x and y, both the 
p (Figure 4.16A) and g (Figure 4.16B) arrows are defined. 
The p perturbation arrows have been made negative with 
mespect to the definition of Schmucker (1970) (that is, 
they now define the anomalous internal current superimposed 
on the normal currents flowing in the + y direction. If 
these two arrows are added vectorially (p + q) as in 
Figure 4.16C, the spatial extent of the anomaly is out- 
lined by arrows polarized at an angle of approximately 45° 
clockwise to the positive x axis. 

In the two-dipole case the induction arrows 
(Figure 4.16D) are no longer polarized in the y direction. 
Over the anomaly the induction arrows tend to point 
radially away from the center of the anomaly. However, 
away from the anomaly the arrows are polarized approximate- 
ly perpendicularly to the boundary and decrease in magni- 
tude near the corners of the anomaly. This is due to the 
high concentration of circulating anomalous current within 


the conductivity anomaly since the conductivity there is 


+ 
greater by a factor of 3.81 x 10 : compared with the 


Surrounding region. Also, the small size of the anomaly 


contributes to the skewing of the arrows. 
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Figure 4.16 


Horizontal two-dipole source with moment 
polarized in the y direction for dipole (1) 
and in the x direction for the dipole (2). 
The dashed arrows are the in-phase arrows 
and the solid arrows are the quadrature- 
phase arrows. 


(A) p perturbation arrows, 

(B) q perturbation arrows, 

(C) p + gq perturbation arrows, 
(D) induction arrows. 


The surface of the conducting region near 
the anomaly is shown and the inner square 
region represents the anomalous structure. 
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For comparison the perturbation and induction 
arrows of two orthogonal uniform sources polarized in the 
x and y directions are presented in Figure 4.17. As found 
in the single dipole case, the p, q and p + q perturbation 
quadrature-phase arrows are greater in magnitude compared 
with the in-phase arrows for the uniform field source. 
Also, as before, the perturbation arrows are more perturbed 
in the uniform field case. The quadrature-phase induction 
arrows are greater in magnitude compared with the in-phase 
arrows for the Wiaoen Source field case. The orientation 
of the induction arrows is very similar. 

As pointed out by Cochrane and Hyndman (1970), 
the problem of defining the normal field for the calcula- 
tion of induction arrows can be partially overcome by 
substituting the total field as the normal field in the 
transfer function calculations, since this 1s easier to 
apply in the analysis of experimental observations. In 
their calculations Cochrane and Hyndman (1970) assumed 
that the magnetic field component in the z direction was 


totally anomalous which thus allows a calculation of the 


induction arrows. This assumption cannot be made when 


calculating the anomalous field for the horizontal 


components of the magnetic field. Figure 4.18 is a 


recalculation of the arrows in Figure 4.16 using the 


Cochrane and Hyndman (1970) suggestion. As can be seen 
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Figure 4.17 As Figure 4.16 but for two orthogonal 
uniform sources with the normal magnetic 
field polarized in the x and y directions. 
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Figure 4.18 


Horizontal two-dipole source with moment 
polarized in the y direction for dipole (1) 
and in the x direction for dipole (2). The 
transfer functions are calculated as 
suggested by Cochrane and Hyndman (1970). 
The dashed arrows are the in-phase arrows 
and the solid arrows are the quadrature- 
phase arrows. 


(A) p perturbation arrows, 

(B) gq perturbation arrows, 

(C) p + gq perturbation arrows, 
(D) induction arrows. 


The surface of the conducting region near the 
anomaly is shown and the inner square region 
represents the anomalous structure. 


232 


ewe 


v 
i 
¢ 
i 
e 
Sissy) ener “ 
: Ha a JTS 4 d | 4 (5 . v fee 


ee ae er a 

a = = 7 
a a < aN a Nae Sie hs ee TC one. « ©. ~ — -, r 
= See OE Cg il oT aire ee canes 


A 
4 
A 


a ad 


ea 


| 


Se 


ase can na tem 
a 
. 
: 
/ 


Sores | 


\ 


ee oe 
Ron. Yoosee To Gare v 


¥-. 


., 


ceed aecoenteneee’ geo’ we 


GES eee 


FOS a aie enema Nea - o " 1 eg ee 4 
>, eee ce, 1. ‘ pote. sidieds hth nomen 
| : f Ee es aes ys ; | y di tek ti on 'tér! di pt 
a eae ae Tee 
pacar mata wor dre palguiate 
| “> Soa sgh ade | apd yndpa 
. es. sat e; it n+ Deas 


a in wwaw p ithe 9 
[ 


{ ' 
ap ——f 4 
ie a} 
=e 
{ 
= 
ft. f 
f ' 
bj 
. - 4 . 
: eee 


| 
| 
t 
i 
t 
t 
a 
aa 
ey waren 
E 
: eae: 


re 

bd 

ie 
Lge 
i 
Lae). 


i 


ad 
y 
' 
' 


——— Ss aes Ss 
DS re 
. 
3 
5 
hd re 
b - <<. 
¢ a r 
Do 
a 
a 
= 
- > 
a 
= 
ag 
* 


2 of the conduetans region W : 
ae and the Taner: RP bates 8! , 
“Av aaTe us 
Oh att et ete Boy * Sees ee 


es oa raps Ye 


a 


me thes 


i eS es Bee. 


ial 
= 


NNN 


240 


in Figure 4.18, little change is noted in the perturbation 
arrows, which are dependent on the horizontal magnetic 
field components. However, the anomalous field used to 
calculate the perturbation and induction arrows in 

Figure 4.18 is exact. 

In this respect, it is of interest to determine 
the effect on the calculations if the anomalous and normal 
fields are approximate. In Figure 4.19 the perturbation 
and induction arrows are given for a calculation in which 
the normal field is assumed to be constant and the anoma- 
lous field is the result of subtracting the normal field 
from the total magnetic field. The vaiue for the normal 
field was obtained by averaging the field values along the 
boundary of the finite difference mesh. The boundary 
forms a square with a perpendicular distance of 60 km from 
the center of the anomaly. This will give a reasonable 
approximation to the normal field since the dipoles are 
located near the anomaly and some symmetry is expected. 
Also the field is reasonably uniform over the region of 
the anomaly. For dipoles placed far from the anomaly a 
better approximation would be a linear interpolation 
across the grid and this latter approximation may also 
be useful for large anomalies. In Figure 4.19 a Slight 


Skewing of the arrows is noted but no appreciable change 


in their configuration is apparent. 
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Figure 4.19 


As Figure 4.18 
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In Figure 4.20 the calculation is done as in 
Figure 4.19 but the total field is substituted for the 


constant normal field. 


4.5 Conclusions 

The transfer functions and related perturbation 
and induction arrows associated with the electromagnetic 
fields induced in various conductivity anomalies by 
uniform, two-dimensional non-uniform and horizontal dipole 
sources oriented in various directions have been studied 
in this chapter. For single sources the induction arrows 
are polarized in a single direction, iach. foe the uniform 
and two-dimensional non-uniform sources, is in the direc- 
tion perpendicular to source current flow. This is due to 
the polarization of the normal field in the uniform source 
case and the polarization of the horizontal normal field 
component and the smallness of the vertical normal field 
component in the two-dimensional non-uniform case. For 
the horizontal dipole the induction arrows are polarized 
in the direction of the dipole moment. In the region 
immediately beneath a horizontal dipole source the normal 
magnetic field component, which is in the direction of the 
dipole moment; is much greater in magnitude than either 
the vertical or perpendicular normal magnetic field 


component, hence the polarization of the induction arrows. 
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Figure 4.20 As Figure 4.18 but for averaged normal field 
values with substitution as suggested by 
Cochrane and Hyndman (1970). 


| 


Seesacc 
~ 


me 


| 
| 


| 


Fenn ean 


ase P oo Foon 


245 


Oe TE | ~<a ee ken SES a ae ae 


cbf 
i 
' 


iimestsa =e an guilds Simlbich wile” eh easy ide 


eee ee - 


{ 
4 hus 9 Yor ‘efetsy 
eefrnanen (1970) F 


ae. et ‘he 
a of LOA ee ee 


\ 
5 jp een \% 


yi * 
yi 
an 
2 
7 ee 
| ' 
' 
ev 
~ a 
4 
a) 


246 


The perturbations arrows for the Single sources for which 
only the q arrow is defined qualitatively resemble the 
anomalous current outside the anomaly. However, near and 
over the anomaly, the in-phase induction arrows increase 
markedly in magnitude regardless of whether the anomaly is 
more conductive or resistive than the surrounding material 
which would not be expected of the anomalous current 
Vectors. 

Various methods were used to calculate the trans- 
fer functions for the multiple source orientation models 
presented. The induction arrows are no longer polarized 
and now point toward (or away from) the conductivity anomaly. 
born p and g perturbations are defined for Chese cases. 

In the uniform source case the concept of a 
normal field-following Schmucker (1970) is used in the 
transfer function calculations and is compared with the 
Cochrane and Hyndman (1970) transfer function calculations. 
It is found that perturbation and induction arrows are 
stable with frequency for the range considered. For this 
source the introduction of the Cochrane and Hyndman (1970) 
substitution produces no appreciable change in the 
configurations of the perturbation and induction arrows. 

It is found that the perturbation arrows when summed 
vectorially do not compare well with the anomalous current 


vectors calculated, nor do they compare with the calculated 
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electric field vectors. However, this new (pq) 
perturbation arrow is very sensitive to changes in 
conductivity and can be used to outline the Spatial 
extent of the anomaly much in the same way as one would 
eserelectric field information 1f it could be accurately 
obtained. The degree of randomness of the source 
orientations was also investigated. Less than five 
source field polarizations can be used if the sources are 
Wagwedst JO" spatially apart. Small angular spreads of 
source field polarization produce skewing of the 
induction arrows. 

Similar calculations were done for the two- 
dimensional non-uniform sources and the horizontal dipole 
sources. In addition to the Schmucker (1970) and Cochrane 
and Hyndman (1970) formulation of the transfer functions, 
the effect of both an approximate anomalous and normal 
field was considered. By using an averaged field value 
for the normal field, the anomalous field is determined 
by subtracting the averaged normal field from the total 
field. The induction arrows are skewed slightly for the 
Cochrane and Hyndman (1970) calculation due to the source 
effect for both source types with the horizontal dipole 
induction arrows showing slightly more deflection due to 
the increased source effect. The averaged normal field 


calculation shows increased skewing for both source types. 
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in general, though,it can’ be’ concluded that the induction 
arrows are stable for a variety of anomalous and normal 
field calculations. The perturbation arrows were more 
stable than the induction arrows for these changes and 
could again be used to outline the spatial extent of the 
be nenarye In the two-dimensional non-uniform source 
calculations the effect of a small angular spread in 
source field polarizations showed the severe skewing 


that“can result. 
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CHAPTER 5 
CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH 


A study of the three-dimensional local induction 
problem has been presented in this work. Solutions for 
the electromagnetic fields were obtained for various 
sources and laterally non-uniform conductivity distribu- 
tions. These calculations were used to compare with 
analogue model measurements for three different source 
configurations. Transfer functions and their related 
perturbation and induction arrows were also calculated for 
the three source types and various embedded conductivity 
anomalies. 

The local induction problem for embedded 
three-dimensional conductivity inhomogeneities for various 
sources was solved through the use of the finite 
difference technique of Lines and Jones (1973). Solutions 
for the three source configurations (uniform, two- 
dimensional, non-uniform and horizontal dipole) above a 
uniform layered conductor were obtained and used as 
boundary conditions for the finite difference mesh. In 
the case of the two-dimensional non-uniform source a 
computer program was developed to calculate the 
electromagnetic fields associated with any arbitrary 


two-dimensional source intensity current distribution 
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over a layered earth. Utilizing the calculated boundary 
conditions, three source field conductivity structure 
combinations were considered. 1. A uniform current source 
"above an eyed structuner in» thea deepy oceans) 2% 4 Ae Wa tichu- 
dinally confined ionospheric electrojet (across which 
current is uniformly distributed) above L-shaped sulfide 
ore body embedded in material with an average continental 
conductivity; 3. A horizontal dipole above an ore body of 
relatively low conductivity which is embedded in poorly 
conducting rock. The profiles of selected electromagnetic 
field quantities and appropriate phases for various 
positions relative to the anomalous structures were 
presented as well as three-dimensional amplitude figures 
of the eee rde and:magnetic field components. The effect 
of the Cae was clearly apparent in all profiles except 
that of the apparent resistivity. The effect of charge 
concentration at the boundary of the anomaly is evident 

in the |E, | and p, component profiles furthest from the 
anomaly. The p, curve exhibits a slight anomaly in the 
opposite sense to the change which would be expected. 

This effect is due to the charge on the boundary and 
affects the p, calculationwiori thisrpeofi le. the: three- 


dimensional amplitude figures of the six electromagnetic 


field components also exhibits strong source effects. 


The corners of the anomalies could easily be identified by 


amplitude maxima in JH I- The sides of the anomaly were 
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marked by increases in the amplitude of [HI. 

A comparison of theoretical calculations with 
analogue model measurements was made for three source 
field configurations over an ore body of relatively low 


conductivity embedded in poorly conducting rock. 


.. Calculations and measurements were made for uniform, 


two-dimensional non-uniform and horizontal dipole sources. 
Selected electromagnetic field quantities and their 
phases were presented in profiles for various positions 
relative to the conductivity structure. It was found that 
in the finite difference method used a transition zone 
must be introduced in the conductivity configuration if 
both high frequency and high conductivity contrast are 
used. Even with the refinement differences were noted in 
the three-dimensional electromagnetic field profiles for 
the theoretical and analogue models over the transition 
zone. In spite of these limitations a good degree of 
compatibility was demonstrated for the two methods and 
their simultaneous use in studying the effects of 
Gahtudedy tiny anomethes is indicated. In addition it was 
found that ene chia ceed’ source calculations wel] 
represented the analogue model sources. 

Transfer functions and related perturbation 
and induction arrows associated with the theoretical 
electromagnetic field calculations were computed. The 


uniform, two-dimensional non-uniform and horizontal 
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dipole sources oriented in various directions and posi- 
tioned variously with respect to the conductivity struc- 
ture were used. A number of methods (i.e., Schmucker 
(1970); Cochrane and Hyndman (1970)) were used to calcu- 
late the transfer functions. Depending on the method 
used and the approximations made for the normal and 
anomalous fields a slight amount of skewing in the arrows 
was noted. In this respect the greatest skewing noted 
was due to the use of averaged normal and approximate 
anomalous fields. However, in general, the induction 
arrows were stable for the variety of anomalous and 
normal field calculations used. The perturbation arrows 
were even more stable in this respect. Little source 
effect was noted for either induction or perturbation 
arrows. Only if a very small angular spread in source 
field polarizations was used would severe skewing be 
noted. This was most evident for the two-dimensional 
non-uniform sources far from the anomaly. As would 
possibly be the case for measurements made at mid- 
latitudes due to auroral electrojets. The results showed 
that the induction arrows could still be used to indicate 
the flow of anomalous currents and point in the direction 
of conductivity anomalies. The perturbation arrows when 
added vectorially (p + q) were not found to compare well 


with either the anomalous current or the anomalous 
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electric field vectors. However, this new arrow was 
found to be very sensitive to conductivity changes in 
much the same manner as the electric field and could be 
used to outline the spatial extent of the anomaly. 

Whe difficuity.of-definingia:nonmal.field, in 
the perturbation and induction arrow transfer function 
analysis emphasizes the desirability of defining the 
transfer functions in terms of the observed surface field 
components. A correlation between the spatial deriva- 
tives of the horizontal and vertical magnetic field 
components exists. Arrows could be defined from the 
calculated transfer functions such that the boundaries 


of the anomalous conductivity region could be defined. 
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APPENDIX A 


THE CALCULATION OF THE ELECTROMAGNETIC FIELDS OF A 
SHEET CURRENT SOURCE WITH ARBITRARY SPATIAL INTENSITY 
DISTRIBUTION OVER A LAYERED HALF SPACE - 

THE COMPUTER PROGRAM AND ITS APPLICATION 


A.1 Description of the Program 


The general two-dimensional program for a sheet 
current source with an arbitrary intensity distribution is 
composed of two programs in which the solution for a general 
Source is constructed from a number of elementary sources 
by superposition. The first program calculates the field 
components over a two-dimensional grid with a layered sub- 
Surface due to a current source with a Gaussian spatial 
intensity distribution of small half width (termed an 
elemental Gaussian). The second program combines the solu- 
tions for the field components of a number of spatially 
shifted elemental Gaussian sources such that a rectangular 
distribution of current intensity is closely approximated. 
With an appropriate change in parameters the second program 
is also used to approximate the arbitrary current intensity 
distribution in a piecewise continuous manner by weighting 
a number of spatially shifted rectangular current distribu- 
tions in accordance with the trapezoidal rule. The solution 


for the field components due to the arbitrary source and the 


layered subsurface is thus obtained. 
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A.2 The Program ELEMENTAL GAUSSIAN 

The ELEMENTAL GAUSSIAN program consists of a main 
program and two subprograms, INTEG and CONST. The main 
program listing is given in Figures A.1 and A.2. The 
subprograms INTEG and CONST are given in Figures A.3, A.4 
and A.5, A.6 respectively. 

The main program calculates the field component 
values at a specified number of points (NT) along a level 
at a constant height for positive y values for a Gaussian 
current distribution at z, = -|h|. 

The integrations for ap H" and Ho are performed 


ay 
through the call to INTEG. The integration subroutine 


 (INTEG) calls the subroutine CONST in which the integrands, 


excluding the terms cos(sy) for ER. and HY and sin(sy) for H, 
(see Chapter 2), are calculated for all values of s at the 
Specified level. For ES the integrand is: 
=0 Zz 
(eh oe: ete Be” }G(si> ss costsy jenms (Pat) 
For Hy the integrand is: 


ae 0,2 
(-i/wy,){8, [Ape ) - Bre © ]e(s)} + cos(sy) 


(A.2) 
For ue the integrand is: 


-9 Z 0.2 
(i/ou,) {ts [Ae n+ Be” Jce(s)} + sin(sy) .  (A.3) 
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¥HSHEAKKKE ELEMENTAL GAUSSIAN*#* tH He RK 


CAUTION CAUTION 


THE PARAMETERS Z21,D1,D2,P9,20,HW ARE EXPRESSED IN METERS 
THE VALUES IN THE ARRAY Z ARE EXPRESSED IN KILOMETERS, 


Z1---HEIGHT OF UPPER BOUNDARY ABOVE SURFACE (A NEGATIVE VALUE) 
Di---THICKNESS OF FIRST LAYER---IF HALF SPACE IS DESIRED SET DT 
GREATER THAN 1.0E+6 METERS, 

D2---THICKNESS OF SECOND LAYER---IF ONLY ONE LAYER IS USED SET D2 
GREATER THAN 1.0E+6 METERS---LIKEWISE FOR A HALF HALF SPACE, 
OI-~--MAXIMNUM INTENSITY OF THE GAUSSIAN CURRNT DISTRIBUTION CENTERED 
AT THE ORIGIN. 

ZO---HEIGHT OF CURRENT SOURCE ABOVE SURFACE (A NEGATIVE VALUE) 
HWe---HALF WIDTH OF GAUSSIAN CURRENT DISTRIBUTION 

Ni---NUMBER OF LEVELS TO BE CALCULATED 

N2,N3,N4,---THE MEANING OF THESE VALUES IS EXPLAINED IN FORMAT 
STATEMENT 112. 

NB---NUMBER OF KILOMETERS BETWEEN CALCULATED FIELD VALUES 
NT~--TOTAL NUMBER OF POINTS PER LEVEL 

NS---ONE HALF THE NUMBER OF INTEGRATION INTERVALS, THIS NUMBER 
MUST BE AN EVEN NUMBER, 


THE MKS SYSTEM OF UNITS IS USED THROUGHOUT THIS PROGRAM. 


Me he Me Me Me ee ae ap ee ae 2 Bee Me ah ee ae ae eae ae tafe af ae he afc i he ie ae ahs ae 


INTEGER Z(25) 

DIMENSION CON (4) ,¥ (2001) 

COMPLEX ODD(2001) , EVEN (2001) ,ORIGIN (2001) , EX (2001) 

COMPLEX ODDHY (2001) ,ODDHZ(2001) , EVENHY (2001) , EVENHZ (2001) , ORIGHY (2 
1001), ORIGHZ (2001) , HY (2001) , HZ (2001) 

COMMON/A1/0I,0M,UM,HW,PI,20, A,B,D1,D2,NT,N5,N8 

READ(5,100) Z1,D1,D2,01,20,HW,N1,N2,N3,N4,NB 

READ(5,101) NT,NS 


Ze----THIS ARRAY CONTAINS THE INTER LEVEL SPACINGS. THE FIRST 
SPACING IS THE DISTANCE BETWEEN THF FIRST LEVEL AND THE TOP BOUNDARY. 


THESE VALUES ARE IN KILOMETERS. 


IF (N1.EQ.1) GO TO 103 
READ(5,102) (Z(K) »K=1,N1) 
GO TO 104 

READ(5,102) Z(1) 


-=--FREQUENCY IN HERTZ. - 
eee OF INTEGRATION. INTEGRAND IS ASSUMED NEGLIGIBLE AT B. A=0. 


READ(5,105) FREQ,B 


HE VARIOUS LAYERS 
Tae CONTAINING THE CONDUCTIVITIES OF T 

Eh peepeen THE AIR LAYER (CONDUCTIVITY ZERO) AND PROCEEDING DOWNWARD. 
FOUR CONDUCTIVITIES MUST BE SPECIFIED. IF A HALF SrA’ eee 
IS DESIRED , THE THREE SUBSURFACE CONDUCTIVITIES MUST BE SPEC 


TO BE THE SAME VALUE. 


Figure A.| ELEMENTAL GAUSSIAN main program. 
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READ(5,106,END=107) CON 
107 A=0.0 

PI=4. O*ATAN(1.0) 

OM=2. 0*PI* FREQ 


UM---PERMEABILITY OF FREE SPACE 


aqanNnAN 


UM=1. 262-6 

WRITE (6,108) 

WRITE (6,109) 21,D1,D2,01,Z0,HW,N1,N8 
WRITE(6,110) FREQ,B 

WRITE(6,111) (CON(I),I=1,4) 

J2=N2=1 

J3=N3-1 

WRITE(6,112) NU,J2,33 


SELECT LEVELS AND DO INTEGRATION FOR FIELD VALUZS BY CALLING INTEG 


aaNAaN 


113 I=1 
Z2=Z1+FLOAT (Z (I) ) *1.0000E+#03 
CALL INTEG (ODD,ODDHY,ODDHZ, EVEN, EVENHY ,EVENHZ, ORIGIN, ORIGHY, ORIGHZ 

1,EX,HY,HZ,Z2,CON,Y) 
WRITE(6,114) I 

WRITE (7) EX 

WRITE(8) HY 

WRITE (9) HZ 

WRITE (6,115) 

WRITE (6,116) (EX (J) ,J=N2,N3,N4) 
WRITE (6,117) 

WRITE (6,116) (HY (J) ,J=N2,N3,N4) 
WRITE (6,118) 

WRITE (6,116) (HZ(J) ,J=N2,N3,N4) 
I=1I+¢1 

IF (Ie LE.N1) GO TO 113 

STOP 

100 FORMAT (6E10.4,12,14,14,14, I4) 

101 FORMAT (I4,15) 

102 FORMAT (213) 

105 FORMAT (P7.5,£10.4) 


106 FORMAT (E10.4 
108 FORMAT (* 1! 21", 9X,4D1", 9X, "D2", 9X, "OL" ,9X,"ZO", 9X, "HW! , 9X, NI, GX 


1,°N8*) 
109 FORMAT ('-",6(E10.4,1X) ,12,2X,I4) 
Mb Sr ecrnns =",F7.5," : UPPER LIMIT OF INTEGRATION =", E10. 4) 


f—9 9 IVITIES= ',4(E10. 4, 2X) ) 
Ws Beer at iconeuet FIELD oe ARE GIVEN AT INTERVALS OF *,14, 1X 
1,°KM. OVER A RANGE OF ',I4," KM. TO si 8 ule) 
114 FORMAT ('-',"LEVEL',I4,° COMPLETED‘) 
115 FORMAT ('-',*EX'! 
116 Be aruach (101.5 (8106 471%, 210.8, 3x19) 
117 FORMAT('-", HY") 
118 FORMAT (*-','HZ') 
END 


Figure A.2 ELEMENTAL GAUSSIAN main program (continued). 
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SUBROUTINE INTEG (ODD, ODDHY ,ODDHZ, EVEN, EVENHY, EVENHZ, ORIGIN, ORIGHY, 
1ORIGHZ,EX,HY,HZ,Z2,CON,Y) 


THIS SUBROUTINE DOES THE INTEGRATION THROUGH THE USE OF FILON*'S METHOD 


DIMENSION CON (4) ,Y (NT) 

COMPLEX ORIGIN (NT) ,ODD (NT) ,EVEN(NT) , EX (NT) 

COMPLEX ORIGHY (NT) ,ORIGHZ (NT) ,EVENHY (NT) , EVENHZ (NT) ,ODDHY (NT) ,ODDH 
1Z(NT) ,HY (NT) ,HZ(NT) 

COMPLEX CONSTT,CONSHY,CONSHZ,ELIN, HYLIN, HZLIN 
COMMON/A1/01,0M,UM,HW,PI,Z0,A,B,D1,D2,NT,N5, NB 


S---INTERVAL OF INTEGRATION, 
S= (B-A) / (FLOAT (N5) *2.) 


Yr~~AFRAY CONTAING THE DISTANCES FFOM THE ORIGIN IN THE POSITIVE Y 
DIRECTION AT WHICH FIELD CALCULATIONS ARE PERFORMED. 


DO 200 I=1,NT 

Y (I) =FLOAT (I-1) *1000. 
ORIGHY (I) =(0.0,0.0) 
ORIGHZ (I) = (0.0,0.0) 
EVENHY (I) =(0.0, 0.0) 
EVENHZ (I) = (0.0,0.0) 
ODDHY (I) =(0.0, 0.0) 
ODDHZ (I) = (0.0, 0.0) 
ORIGIN (I) =(0.0,0.0) 
EVEN (I) = (0.0,0.0) 
ODD(I)=(0.C, 0.0) 


ARRAYS PREFIXED WITH ODD CONTAIN THE CALCULATIONS OF THE INTEGRAND 
PERFORMED FOR THE ODD NUMBERED INTERVALS OF INTEGRATION, IN ACCORDANCE 
WITH FILON'S METHOD. 


ARRAYS PRFFIXED WITH EVEN CONTAIN THE CALCULATIONS OF THE INTEGRAND 
PERFORMED FOR THE EVEN NUMBERED INTERVALS OF INTEGRATION, IN ACCORDANCE 
WITH FILON'S METHOD. 


T=S 

DO 201 J=1,N5 

CALL CONST (T,22,CONSTT,CONSHY,CONSHZ,CON) 
DO 202 I=1,NT, N8 

ODDHY (1) =ODDHY (I) ¢CONSHY*COS (T*Y¥ (T) ) 
ODDHZ (I) =OLDHZ (I) +CONSHZ*SIN (T*Y (I) ) 

ODD (I) =ODD (I) +CONSTT*COS (T*Y (TZ) ) 

T=T+S¢S 

NB=N5=-1 

G=S+S 

DO 203 J=1,NB 

CALL CONST (G, 22, CONSTT, CON SHY ,CONSHZ,CON) 
DO 204 I=1,NT,N8 

EV ENHY (I) =EVENHY (I) *CONSHY*COS (G*¥ (T)) 
EVENHZ (I) =EVENHZ (I) *CONSHZ*S IN (G*¥ (T)) 
EVEN (I) =EVEN (I) *CONSTT*COS (G*Y (1) ) 


Figure A.3 Subroutine INTEG. 
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G=G#S¢S 

V=1.000E-20 

CALL CONST (V,22,CONSTT,CONSHY,CONSHZ,CON) 
DO 205 I=1,NT, NB 

ORIGHY (I) =CONSHY 

ORIGHZ (I) =CONSHZ 

ORIGIN(I) =CONSTT 


FILON*S INTEGRATION METHOD 


THE=Y (I) *S 

IF(THE .LE. .1) GO TO 206 

ALP= (THE**2¢THE*SIN (THE) *COS (THE) = (2.* (SIN (THE) ) **2)) /THE 
BETA=2, * (THE* (1. ¢( (COS (THE) ) **2) )-2. *SIN (THE) *COS (THE) ) /THE**3 
GAMMA=4, * (SIN(THE) ~THE*COS (THE) ) /THE**3 

GO TO 207 

BETA=2./30 +2. * (THE**2) /15.-4.,* (THE**4) /105,42.* (THE**6) /56 76 
GAMMA=4./3.-2.*(THE**2) /15. 4 (THE®*4) /210.> (THE**6) /11340. 
ALP=( (2. *THE**3) /45.) = ((2e*THE®*5) /315.) +( (2. *THE*®*7) /4725.) 
HY (I) =S* (BETA* (EVENHY (I) -. 5*ORIGHY (I) ) *+GAMMA*ODDHY (I) ) 

HZ (I) =S* (ALP*ORIGHZ (I) +BETA*EVENHZ (I) +GAMMA*ODDHZ (I) ) 

EX (I) =S* (BETA* (EVEN (I) =. 5*ORIGIN (I) ) ¢+GAMMA*ODD (I) ) 


LINEAR INTERPOLATION BETWEEN CALCULATED POINTS IN A GIVEN LEVEL. 


N9=N8+1 

N10=N8-1 

DO 208 I=N9,NT,NB 
ELIN=(EX (I) -EX (I-N@) ) /FLOAT(N8) 

HY LIN= (HY (I) “HY (I-N8) ) /FLOAT (N8) 
HZLIN= (HZ (I) -HZ(I-N8) ) /FLOAT (N@) 
DO 208 J=1,N10 

EX (I-N8¢J) =EX (I-N8) ¢ELIN* FLOAT (J) 
HY (I-N8¢J) =HY (I-N8) *HYLIN*FLOAT (J) 
HZ (I-N8+J) =HZ (I-N8) *HZLIN* FLOAT (J) 
RETURN 

END 


Figure A.4 Subroutine INTEG (continued). 
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SUBROUTINE CONST (S,22,E1,HY1,HZ1, CON) 
DIMENSION CON (4) 
COMPLEX TH (3) -W(3) ,Q, CEXP,R, B1,CSQRT,A2,A3,A4,B2,B3,01,21 
COMPLEX HY1,HZ1 
COMPLEX CPRIM,CMPLX 
COMMON/A1/01,0N,UM,HW,PI,20,A,B,D1,D2,NT,N5, NB 


THIS SUBROUTINE CALCULATES THE INTEGRAND VALUES, 
THE SUBSURFACE CAN CONTAIN A MAXIMUM OF TWO LAYERS AND ONE HALF SPACE, 


CPRIM=(0.0,-1.0) * ( (OL*ON*UM*HW) /(SQRT(2.0*PI) *S) ) *EXP(S*Z0-S** 2* HW 
1*42/2.) 
Z=S**2 
DO 300 N=1,3 
Y1=O8*UM*CON (N+1) 
W(N) =CMPLX (Z, Y1) 
TH (N) =CSQRT(W(N) ) 
IF (D1.GT.1.0E+6) GO TO 301 
IF (D2.GT.1.0E+6) Q1=CMPLX(0.0,0.0) 
IF (D2.GT.1.0E+6) GO TO 301 
Q1= ((TH (2) -TH (3) ) / (TH (2) +TH (3) )) *CEXP (-2.*TH (2) *D2) 
IF (D1.GT.1.0E+6) GO TO 302 
Q= ((TH(1) #TH (2) ¢(TH(1) -TH(2) ) *Q7*CEXP (2. *TH(2) *D1)) / (TH (1) -TH (2) + ( 
1TH (1) +TH (2) ) *Q7*CEXP(2.*TH (2) *D1))) *CEXP (2.*TH(1) *D1) 
R= (S* (Q¢1.)) / (TH (1) * (Q-14)) 
BI=(R-1.)/(R+1-) 
IF (D1.GT.1.0E+6) B1=(S=-TH(1))/(S¢#TH (1) ) 
IF(Z2 .LE. 0.0) GO TO 304 
IF (D1.GT.1.0E+6) AY=CMPLKX(1.0,0.0) +B1 
IF (D1.GT.1-0E+6) GO TO 307 
A2=((1.+B1) *Q) /(Q*1.-) 
B2=A2/0 
IF(Z2 -LE. D1) GO TO 305 
A3=( (A2*+B2*CEXP (2. *TH (1) *D1) )/ (1. +Q 1*CEXP (2. *TH (2) *D1) ) ) *CEXP( (TH ( 
12) -TH (1) ) *D1) 
B3=A3¥*Q1 
IF (D2 .GT. 1.0 E#6) GO TO 303 
IP(Z2 .LE. D2) GO TO 306 
A4= (A3+B3*CEXP (2*TH (2) *D2) ) *CEXP ( (TH (3) ~TH (2) ) *D2) 
GO TO 307 
A4=A3 
GO TO 307 Peg teri 
= -S*Z2) +B1* EXP (S* 
HYIe((OL0 -1.0) 7 (ON#UM) ) #5 (EXP (-S#Z2) “BIEXP (S#Z2) ) #CPRIM 
HZ1= (((0.0,1.0)*S) / (OM*UM) ) *E1 
GO TO 308 
= * -TH (1) *Z2) ¢B2*CEXP (TH (1) *Z2) ) *CPRIM 
Faroe 4 0) ) lomennyy *R (1) *(A2*CEXP (<TH (1) *Z2) -B2*CEXP (TH (1) *Z 
12) ) *CPRIM 
HZ1=(( (020,120) *S) /(OM*UM) ) *E1 
GO TO 308 
= * - 2) *Z2) +B3*CEXP (TH (2) *Z2) ) *CPRIM 
HY Te (Oc O oe 400) 2 (OMMUM) ) #TH (2) # (ABHCEXD (—TH (2) #22) ~B3#CEXP (TH (2) #2 


12) ) *CPRIM 


Figure A.5 Subroutine CONST. 
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HZ1=(((0.0,160) *S) / (OM*UM) ) ¥E1 
GO TO 308 


307 E1= (AY*CEXP (-TH (3) *Z2)) *CPRIM 


308 


HY 1=((0.0,-160) /(OM*UM) ) *TH (3) * (AG*CEXP(-TH (3) *Z2)) *CPRIM 
HZ 1=(( (020,120) *S) /(OM*UM) ) *E1 

RETURN 

END 


Figure A.6 Subroutine CONST (continued). 
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C(s) is as given in equation(2.94) of Section 2.5.Upon 
returning to INTEG the integration is performed by the use of 
Filon's method (Tranter, 1965). The main program then writes 
the electric and magnetic field values for that level on a 
file or magnetic tape. This process is then repeated for 


Subsequent levels. 


A.3 The Program RECTAZOID 

This program sums a finite number of elemental 
Gaussian solutions to obtain an approximate solution for a 
rectangular current distribution. The user specifies the 
number of elemental Gaussians to approximate the rectangular 
Current distribution and assigns their source coefficients 
the magnitude one. The elemental Gaussian solution for 
positive values of y calculated by the previous program is 
read into both positive and negative y regions. In this, the 
vertical magnetic field quantities must be made negative in 
the negative y region since the vertical magnetic field is 
asymmetric about the origin. This solution is then shifted 
and summed to approximate the rectangular current solution. 
The result is then stored. This program is again used to 
give the solution for a piecewise continuous source composed 
of a summation of rectangular current sources. The general 
Source configuration must digitized by the user at regular 
intervals (equal to the width of the rectangles) and these 


Source coefficients are then supplied to the program in the 
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array SA. The general solution is obtained by superimposing 
rectangular source solutions which are shifted with respect 
to one another and multiplied by the appropriate source 
coefficients in accordance with the trapezoidal rule. The 
final solution is thus obtained for the layered earth and 

is that of a piecewise continuous source which approximates 
a continuous source. The listing of RECTAZOID is given in 


Figures A.7 and A.8. 


A.4 Various Parameters and Their Pigeers on the Solution 
In the programs many parameters can be varied and 

the values of these can affect the accuracy of the final 

solution. It is important to consider these parameters and 


the reasons for choosing their values. 


A.4.1 Elemental Gaussian 


The final result depends greatly upon the accuracy 
of the elemental Gaussian solution. The important 


parameters here are B, N5 and N8. 


For half widths of approximately 10 km the source 
SZ 


fe) 
function C(s) is governed by the terme ~. Therefore, B 
should be taken to be at least ten times the damping constant 
my} z |) for a source in free space (see Hermance and Peltier 
0) 


-4 
(1970). For a layered-earth model a value of 0.8 x 10 © or 


larger will be satisfactory. 
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HOS eEHEAKRECTAZOL D** #4 KM 


NT---NUMNBER OF FIELD POINTS USED IN THE ELEMENTAL GAUSSIAN PROGAM. 
NLEVEL~--NUMBER OF LEVELS. 


NSOUR---NUMBER OF SOURCE ELEMENTS. MUST BE ODD NUMBER FOR RECTANGULAR 
SOURCE. 

NSC---NUMBER OF THE SOURCE ELEMENT WHICH IS POSITIONED OVER THE ORIGIN. 
NSHIFT---NUMBER OF ARRARY ELEMENTS THAT FIELD IS TO BE SHIFTED FOR EACH 
CURRENT ELEMENT. 

NHZ---IF HZ MAGNETIC FIELD DATA IS BEING RUN FOR A RECTANGULAR SOURCE 
NHZ=1. AT ALL OTHER TINFS NHZ=0. THIS IS TO ALLOW FOR NON-SYMMETRY OP HZ 
MAGNETIC FIELD VALUES. 

THE DIMENSION OF THE SOURCE ARRAY (SA) MUST EQUAL NSOUR. 

NHW---HALF WIDTH OF ELEMENTAL GAUSSIAN. 


We Hee whe ee he he eee ee aie ae ae ae hee eee ae ea Ok 


COMPLEX ET(4001) ,EX(2001) ,EA (4001) 

DIMENSION SA(14) 

READ(5,100) NT,NLEVEL,NSOUR, NSC, NSHIPT,NHZ,NHW 
READ(5,101) N4,N2,N3 

WRITE (6,102) 

WRITE(6,103) NT,NLEVEL,NSOUR,NSC,NSHIFT,NHZ, NHW 
J2=N2+¢1 

J3=N3-1 

WRITE (6,104) N4,J2,J3 

NT1=NT=1 

NT2=NT+1 

NTOT=2*NT=1 

IF (NSHIFT.GT.NHW) GO TO 106 


SOURCE ARRAY FOR RECTANGULAR SOURCE. 


DO 105 I=1,NSOUR 
SA (I) =1.000 
GO TO 107 


SOURCE ARRAY FOR LINEARLY APPROXIMATED SOURCE. 


READ(7) SA 
11=1 
IF (NSHIFT.GT.NHW) GO TO 112 


ELEMENTAL GAUSSIAN FIELD VALUES 


READ (8) EX 


FORM SYMMETRIC FIELD ARRAY (EXCEPT IN CASE OF HZ FOR RECTANGULAR SOURCE) 
FROM ELEMENTAL GAUSSIAN DATA. 


DO 109 I=1,NT 
EA (NT 141) =EX (I) 
EA (I) =EX (NT2-T) 
IF (NHZ.EQ.0) GO TO 111 


Figure A.7 The program RECTAZOID. 
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DO 110 I=1,NT 
110 EA(I) =-EA (1) 
111 IF(NSHIFT.EQ.NHW) GO TO 113 


RECTAZOID FIELD VALUES, 
112 READ(8) EA 


MULTIPLY BY AMPLITUDE OF RECTANGULAR SOURCE CENTERED OVER THE ORIGIN. 


113 DO 114 I=1,NTOT 

114 ET (I) =EA (I) *SA (NSC) 
J=(NSC-1) *NSHIFT 
J1=1 

115 IF (JeZ0.0) GO TO 117 


SHIFT AND ADD FOR SOURCE ELEMENTS IN THE NEGATIVE Y REGION. 


DO 116 I=1,NTOT 
IF ((I¢J).LT.1.OR. (I+J).GT.NTOT) GO TO 116 
ET (1) =ET (I) ¢EA (I¢J) *SA (J1) 
116 CONTINUE 
J=J-NSHIFT 
J1=J31+1 
GO TO 115 
117 K=NSC#1 
J=NSHIFT 
118 IF(K.GT.NSOUR) GO TO 120 


SHIFT AND ADD FOR SOURCE ELEMENTS IN THE POSITIVE Y REGION, 


DO 119 I=1,NTOT 
IF ((I-J).LT.1.OR. (I-J).GT.NTOT) GO TO 119 
ET (I) =ET(I) #EA (I-J) *SA (K) 
119 CONTINUE 
J=J¢NSHIFT 
K=K+1 
GO TO 118 
120 WRITE(9) ET 
WRITE(6,121) 11 
WRITE (6,122) (ET (J) ,J=1,NTOT, NG) 
I1=11+1 
IPF (I1.GT.NIEVEL) GC TO 123 
GO TO 108 
123 STOP 
100 FORMAT (I4, 613) 


101 FORMAT (315 
102 rodnimsee NT NLEVEL NSOUR NSC NSHIFT NHZ NHW*) 


t-1 3,4X, 13 ,2X,13¢5X~¢13, 1X, 13, 1X-13) 
iu Oey ace LepNPLEX FIELD. WALUES ARE GIVEN AT INTERVALS OF *,I5, 1X 
1,°KM. OVER A RANGE OF ',15,* KM. TO ',I5,' KM. ') 
121 PORMAT(*-',*LEVEL',14," COMPLETED‘) 
122 FORMAT (100 (0%,5 (E1004,1X,E10.4, 3X) /)) 
END 


Figure A.8 The program RECTAZOID (continued). 
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A total of 4000 integration intervals (2*N5) is 
sufficient for a half width (HW) of 1-10 km. This value of 
NS improves the accuracy of the apparent resistivity curves 


calculated from 


(Cagniard, 1953) since when calculating the apparent 
resistivity any error in the field values produces a greater 
error in the resistivity values. Satisfactory field values 
may be obtained by using approximately one-sixth this 

value of N5. These values of B and N5 were used in 
Calculating the elemental Gaussian results of Tables 1 and 

2 of Chapter 2. In calculating the 240 km half width 
Gaussian the value of B was taken as 0.1 x 10°. 

The solution for the field values at any particu- 
lar level is calculated at a certain number of points and 
this is indirectly determined by N8. The values for inter- 
mediate points are determined by a linear interpolation 
between the calculated points. The quantity (N8-1) 
represents the number of points between calculated values. 
In all cases except the sample run N8& was set equal to 10. 
If this value is made too large the solution will be less 


accurate, especially when the field solutions are shifted 


and added. 
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‘“A.4.2 RECTAZOID 

The parameters in RECTAZOID which most affect the 
results are NSOUR and NSHIFT. When Synthesizing the 
rectangular source the best approximation is obtained when 
the adjacent elemental Gaussians are one half width apart. 
Because of this, NSHIFT is equal to the half width in 
kilometres. When arbitrary source currents are approximated 
NSHIFT is equal to the width of the rectangles. The width 
of the rectangular source specified by NSOUR can also affect 
the accuracy of the solution. The normalized values of H, 
are most affected by the source geometry and coarse rectan- 
gular approximations to the arbitrary source configuration 


Can noticeably affect these profiles. 


A.5 Sample Run 
Sample calculations have been made and the results 
of this run are given in Figures A.9, A.10, A.11 and A.12. 
The sample run has been designed to provide an 
easy verification of the program. The total running time 
to obtain the results given here is approximately 40 s on 
the IBM 360/67 at the University of Alberta. In the 
ELEMENTAL GAUSSIAN PROGRAM for this example only three 


points are calculated and the remaining ones are obtained 


through linear interpolation. A single value of 50 must be 


Supplied for the Z array since the fields are obtained for 
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the surface level. The value of NT was 2001 and N5 was taken 
to be 20. The other parameters are as specified in Figure 
A.9. Figure A.10 gives the parameters used in RECTAZOID to 
construct the 30-km rectangular source. Since only vertical 
magnetic field data are produced in this run NHZ must be set 
equal to 1 to ensure the non-symmetry of the He field. The 
value of NHZ must be taken as 1 in construction of the 
rectangular source, but when synthesizing the arbitrary 
source from these rectangles NHZ must equal zero. The 
arbitrary source consists of fourteen rectangular sources 
weighted as indicated in the source array SA given in Figure 
A.11. The third rectangular current source is positioned at 
the origin as indicated by NSC. The parameters of Figure 
A.12 produce the vertical magnetic field associated with this 
Source and these field values are given as well in Figure 
A.12. If the electric field or the horizontal magnetic 

field is desired, the elemental Gaussian data for those 
fields could be operated upon by RECTAZOID to produce the 
desired field. In these latter two cases all parameters for 
RECTAZOID would remain the same as before, except for NHZ 


which would be set equal to zero. 
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